feedback calculation
This commit is contained in:
church89
parent
e6e2408e13
commit
75bb200d76
1 changed files with 187 additions and 57 deletions
244
msre.py
244
msre.py
|
|
@ -3,61 +3,114 @@ import openmc
|
|||
import openmc.deplete
|
||||
import openmc.lib
|
||||
import numpy as np
|
||||
from math import log10
|
||||
from math import log10, sqrt
|
||||
from collections import OrderedDict
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
from scipy import stats
|
||||
import math
|
||||
from scipy.optimize import curve_fit
|
||||
import numpy as np
|
||||
import re
|
||||
|
||||
def define_materials():
|
||||
# fuel salt
|
||||
salt = openmc.Material(name="salt", temperature = 911)
|
||||
salt.add_nuclide('Li6',1.31480070E-05)
|
||||
salt.add_nuclide('Li7', 0.262960140146177)
|
||||
salt.add_nuclide('Be9',1.1863E-01)
|
||||
salt.add_nuclide('Zr90',1.0543E-02)
|
||||
salt.add_nuclide('Zr91',2.2991E-03)
|
||||
salt.add_nuclide('Zr92',3.5142E-03)
|
||||
salt.add_nuclide('Zr94',3.5613E-03)
|
||||
salt.add_nuclide('Zr96',5.7375E-04)
|
||||
salt.add_nuclide('Hf174',8.3786E-10)
|
||||
salt.add_nuclide('Hf176',2.7545E-08)
|
||||
salt.add_nuclide('Hf177',9.7401E-08)
|
||||
salt.add_nuclide('Hf178',1.4285E-07)
|
||||
salt.add_nuclide('Hf179',7.1323E-08)
|
||||
salt.add_nuclide('Hf180',1.8370E-07)
|
||||
salt.add_nuclide('U234',1.034276246E-05)
|
||||
salt.add_nuclide('U235',1.009695816E-03)
|
||||
salt.add_nuclide('U236',4.227809892E-06)
|
||||
salt.add_nuclide('U238',2.168267822E-03)
|
||||
salt.add_nuclide('Fe54',2.8551E-06)
|
||||
salt.add_nuclide('Fe56',4.4818E-05)
|
||||
salt.add_nuclide('Fe57',1.0350E-06)
|
||||
salt.add_nuclide('Fe58',1.3775E-07)
|
||||
salt.add_nuclide('Cr50',2.1224E-06)
|
||||
salt.add_nuclide('Cr52',4.0928E-05)
|
||||
salt.add_nuclide('Cr53',4.6409E-06)
|
||||
salt.add_nuclide('Cr54',1.1552E-06)
|
||||
salt.add_nuclide('Ni58',5.8597E-06)
|
||||
salt.add_nuclide('Ni60',2.2571E-06)
|
||||
salt.add_nuclide('Ni61',9.8117E-08)
|
||||
salt.add_nuclide('Ni62',3.1284E-07)
|
||||
salt.add_nuclide('Ni64',7.9671E-08)
|
||||
salt.add_nuclide('O16',5.1437E-04)
|
||||
salt.add_nuclide('O17',1.8927E-07)
|
||||
salt.add_nuclide('O18',9.6440E-07)
|
||||
salt.add_nuclide('F19',5.9409E-01)
|
||||
salt.set_density('g/cm3', 2.3275)
|
||||
def calc_density(temp_c):
|
||||
temp_f = temp_c*9/5 + 32 #Convert from Celsius to F
|
||||
temp_ref = 1181 # F or 638.3 C
|
||||
rho_ref = 2.3275 # g/cm3 @ 1181 F
|
||||
exp_c = -1.18e-4 # 1/F fuel salt expansion coefficient
|
||||
rho_calc = rho_ref * (1 + exp_c*(temp_f-temp_ref)) #adjusted density
|
||||
return rho_calc
|
||||
|
||||
def define_materials(fuel_temp, u235_load=None):
|
||||
# CRITICALITY BENCHMARK, definition based on U235 mass fraction of 0.01408 wt%
|
||||
# fuel salt temperature of 638.3 C and salt density of 2.3275
|
||||
if u235_load is None:
|
||||
rho = calc_density(fuel_temp)
|
||||
tot_mass = 4340.43 * 1000 # Total fuel salt in the loop charge, fixed [g]
|
||||
volume = tot_mass/rho # assuming total volume remains constant
|
||||
|
||||
salt = openmc.Material(name="salt", temperature = fuel_temp + 273.15)
|
||||
salt.add_nuclide('Li6',1.31480070E-05)
|
||||
salt.add_nuclide('Li7', 0.262960140146177)
|
||||
salt.add_nuclide('Be9',1.1863E-01)
|
||||
salt.add_nuclide('Zr90',1.0543E-02)
|
||||
salt.add_nuclide('Zr91',2.2991E-03)
|
||||
salt.add_nuclide('Zr92',3.5142E-03)
|
||||
salt.add_nuclide('Zr94',3.5613E-03)
|
||||
salt.add_nuclide('Zr96',5.7375E-04)
|
||||
salt.add_nuclide('Hf174',8.3786E-10)
|
||||
salt.add_nuclide('Hf176',2.7545E-08)
|
||||
salt.add_nuclide('Hf177',9.7401E-08)
|
||||
salt.add_nuclide('Hf178',1.4285E-07)
|
||||
salt.add_nuclide('Hf179',7.1323E-08)
|
||||
salt.add_nuclide('Hf180',1.8370E-07)
|
||||
salt.add_nuclide('U234',1.034276246E-05)
|
||||
salt.add_nuclide('U235',1.009695816E-03)
|
||||
salt.add_nuclide('U236',4.227809892E-06)
|
||||
salt.add_nuclide('U238',2.168267822E-03)
|
||||
salt.add_nuclide('Fe54',2.8551E-06)
|
||||
salt.add_nuclide('Fe56',4.4818E-05)
|
||||
salt.add_nuclide('Fe57',1.0350E-06)
|
||||
salt.add_nuclide('Fe58',1.3775E-07)
|
||||
salt.add_nuclide('Cr50',2.1224E-06)
|
||||
salt.add_nuclide('Cr52',4.0928E-05)
|
||||
salt.add_nuclide('Cr53',4.6409E-06)
|
||||
salt.add_nuclide('Cr54',1.1552E-06)
|
||||
salt.add_nuclide('Ni58',5.8597E-06)
|
||||
salt.add_nuclide('Ni60',2.2571E-06)
|
||||
salt.add_nuclide('Ni61',9.8117E-08)
|
||||
salt.add_nuclide('Ni62',3.1284E-07)
|
||||
salt.add_nuclide('Ni64',7.9671E-08)
|
||||
salt.add_nuclide('O16',5.1437E-04)
|
||||
salt.add_nuclide('O17',1.8927E-07)
|
||||
salt.add_nuclide('O18',9.6440E-07)
|
||||
salt.add_nuclide('F19',5.9409E-01)
|
||||
salt.set_density('g/cm3', rho_calc)
|
||||
|
||||
else:
|
||||
# REACTIVITY BENCHMARK
|
||||
rho = calc_density(fuel_temp)
|
||||
# Salt composition in kg with 65.25 kg U235 at 649 C
|
||||
salt_comp = {'Li6':507.27*0.00005,'Li7':507.27*0.99995,'Be':293.96,
|
||||
'Zr':513.97,'Hf':0.0029,'U234':0.67,'U235':65.25,'U236':0.27,
|
||||
'U238':141.91,'Fe':0.75,'Cr':0.13,'Ni':0.14,'O':2.27,'F':3103.22}
|
||||
tot_mass = sum(salt_comp.values())*1000
|
||||
# Calculate total volume assuming fuel salt composition at first load reamins fixed
|
||||
volume = tot_mass/rho
|
||||
|
||||
# HEU additions mass fraction
|
||||
heu = {'Li6':4.93*0.00005,'Li7':4.93*0.99995,'U234':0.59,
|
||||
'U235':57.7,'U236':0.24,'U238':3.54,'F':33}
|
||||
|
||||
# Add u235
|
||||
mass_add = u235_load - salt_comp['U235']
|
||||
for comp,mass in salt_comp.items():
|
||||
if comp in heu:
|
||||
salt_comp[comp] += mass_add/heu['U235']*heu[comp]
|
||||
# Recalculate total mass
|
||||
tot_mass = sum(salt_comp.values())*1000
|
||||
|
||||
salt = openmc.Material(name="salt", temperature = fuel_temp + 273.15)
|
||||
for comp, mass in salt_comp.items():
|
||||
splt = re.split(r'\d+', comp)
|
||||
if len(splt) == 2:
|
||||
salt.add_nuclide(comp, mass/tot_mass, 'wo')
|
||||
else:
|
||||
for nuc,frac in openmc.data.isotopes(comp):
|
||||
salt.add_nuclide(nuc, mass*frac/tot_mass, 'wo')
|
||||
# Recalculate density as new total mass and fixed volume
|
||||
salt.set_density('g/cm3',tot_mass/volume)
|
||||
#Get U235 mass
|
||||
salt.volume = volume #Assign volume as Loop total volume in cm3
|
||||
|
||||
#moderator blocks170
|
||||
graphite = openmc.Material(name='graphite',temperature=911)
|
||||
graphite = openmc.Material(name='graphite',temperature=fuel_temp + 273.15)
|
||||
graphite.set_density('g/cm3',1.86)
|
||||
graphite.add_nuclide('C12',1)
|
||||
graphite.add_s_alpha_beta('c_Graphite')
|
||||
|
||||
#inor-8
|
||||
inor = openmc.Material(name='inor-8',temperature=911)
|
||||
inor = openmc.Material(name='inor-8',temperature=fuel_temp + 273.15)
|
||||
inor.set_density('g/cm3',8.7745)
|
||||
inor.add_element('Ni',(66+71)/2,'wo')
|
||||
inor.add_element('Mo',(15+18)/2,'wo')
|
||||
|
|
@ -82,7 +135,7 @@ def define_materials():
|
|||
|
||||
#Control rods inconel clad
|
||||
trace = 0.01
|
||||
inconel = openmc.Material(name='inconel', temperature = 911)
|
||||
inconel = openmc.Material(name='inconel', temperature = 65.6 + 273.15)
|
||||
inconel.add_element('Ni',78.5,percent_type='wo')
|
||||
inconel.add_element('Cr',14.0,percent_type='wo')
|
||||
inconel.add_element('Fe',6.5,percent_type='wo')
|
||||
|
|
@ -110,7 +163,7 @@ def define_materials():
|
|||
inconel.set_density('g/cm3',8.5)
|
||||
|
||||
# SS 316 control rod flexible hose
|
||||
ss316 = openmc.Material(name='ss316')
|
||||
ss316 = openmc.Material(name='ss316', temperature = 65.6 + 273.15)
|
||||
ss316.add_element('C',0.026,'wo')
|
||||
ss316.add_element('Si',0.37,'wo')
|
||||
ss316.add_element('Mn',0.16,'wo')
|
||||
|
|
@ -135,6 +188,7 @@ def define_materials():
|
|||
Al2O3.set_density('g/cm3',3.95)
|
||||
bush = openmc.Material.mix_materials([Gd2O3,Al2O3],[0.7,0.3],'wo')
|
||||
bush.name='bush'
|
||||
bush.temperature = 65.6 +273.15
|
||||
|
||||
#Concrete block
|
||||
concrete = openmc.Material(name='concrete')
|
||||
|
|
@ -192,7 +246,7 @@ def define_materials():
|
|||
mats = openmc.Materials([salt,graphite,inor,helium,inconel,shield,concrete,steel,ss316,sandwater,insulation,bush])
|
||||
return mats
|
||||
|
||||
def build(make_tally=True, plot_geom=True):
|
||||
def build(make_tally=True, plot_geom=True, u235_load=None, fuel_temp=638.3, cr1_pos=46.6, cr2_pos=51, cr3_pos=51):
|
||||
#Clean-up
|
||||
os.system("rm *.xml *.h5 *.out")
|
||||
|
||||
|
|
@ -201,7 +255,7 @@ def build(make_tally=True, plot_geom=True):
|
|||
control_rod1_h5m = 'h5m/msre_control_rod_1e-2.h5m'
|
||||
|
||||
#Materials
|
||||
mats=define_materials()
|
||||
mats=define_materials(fuel_temp,u235_load)
|
||||
|
||||
#Geometry
|
||||
core = openmc.DAGMCUniverse(filename=core_h5m, auto_geom_ids=True, universe_id=1)
|
||||
|
|
@ -216,18 +270,19 @@ def build(make_tally=True, plot_geom=True):
|
|||
cr1_cell = openmc.Cell(name='CR1', region=cr1_region, fill=control_rod1)
|
||||
cr2_cell = openmc.Cell(name='CR2', region=cr2_region, fill=control_rod1)
|
||||
cr3_cell = openmc.Cell(name='CR3', region=cr3_region, fill=control_rod1)
|
||||
#translate control rod1 130 cm (51 inch) up in z direction
|
||||
setattr(cr1_cell, 'translation', [0,0,130+19.2])
|
||||
setattr(cr2_cell, 'translation', [-10.163255,0,130+19.2])
|
||||
setattr(cr3_cell, 'translation', [-10.163255,10.163255,130+19.2])
|
||||
#translate control rods
|
||||
start_pos = 19.2
|
||||
setattr(cr1_cell, 'translation', [0, 0, start_pos + cr1_pos*2.54])
|
||||
setattr(cr2_cell, 'translation', [-10.163255, 0, start_pos + cr2_pos*2.54])
|
||||
setattr(cr3_cell, 'translation', [-10.163255, 10.163255, start_pos + cr3_pos*2.54])
|
||||
geometry = openmc.Geometry([core_cell,cr1_cell,cr2_cell,cr3_cell])
|
||||
|
||||
# Settings
|
||||
settings = openmc.Settings()
|
||||
settings.temperature = {'method':'interpolation','range':(293.15,923.15)}
|
||||
settings.batches = 50
|
||||
settings.inactive = 10
|
||||
settings.particles = 50000
|
||||
settings.batches = 300
|
||||
settings.inactive = 50
|
||||
settings.particles = 100000
|
||||
settings.photon_transport = False
|
||||
source_area = openmc.stats.Box([-100., -100., 0.],[ 100., 100., 200.],only_fissionable = True)
|
||||
settings.source = openmc.Source(space=source_area)
|
||||
|
|
@ -404,7 +459,7 @@ def control_rod_worth(model):
|
|||
pos *= 2.54
|
||||
dx = 4 * 2.54
|
||||
for shift in [-dx,dx]:
|
||||
setattr(cell, 'translation', [0, 0, 19.2+pos+dx])
|
||||
setattr(cell, 'translation', [0, 0, 19.2+pos+shift])
|
||||
res=model.run()
|
||||
with openmc.StatePoint(res) as sp:
|
||||
keff=sp.keff.n
|
||||
|
|
@ -422,6 +477,74 @@ def control_rod_worth(model):
|
|||
plt.ylabel('Reactivity worth')
|
||||
plt.savefig('reac_rod_worth')
|
||||
|
||||
def func (x,a,b):
|
||||
return a*x +b
|
||||
|
||||
def feedback_isothermal(temperatures,u235_load):
|
||||
# temperature arguments in Celsius
|
||||
y_data = []
|
||||
y_err = []
|
||||
for temp in temperatures:
|
||||
model = build(make_tally=False, plot_geom=False, fuel_temp = temp, u235_load=u235_load)
|
||||
res=model.run()
|
||||
with openmc.StatePoint(res) as sp:
|
||||
k=sp.keff.n
|
||||
y_data.append((k-1)/k)
|
||||
y_err.append(sp.keff.s)
|
||||
|
||||
x_data = np.array(temperatures)
|
||||
y_data = np.array(y_data)
|
||||
y_err = np.array(y_err)
|
||||
popt, pcov = curve_fit(func, x_data, y_data)
|
||||
|
||||
plt.figure()
|
||||
plt.errorbar(x_data, y_data, yerr=y_err, label='Sim data')
|
||||
symb = r'$^\circ$C'
|
||||
#propagated error
|
||||
err = r'$\pm${:.4f}'.format(sqrt(sum(y_err**2))*1e5)
|
||||
plt.plot(x_data, func(x_data, *popt), '--', label='{:.3f} [pcm/{}]'.format(popt[0]*1e5,symb))
|
||||
plt.legend()
|
||||
plt.xlabel('Temperature '+symb)
|
||||
plt.ylabel(r'$\rho\,\pm\sigma$')
|
||||
plt.title(r'Isothermal temperature coefficient, {u235_load} kg $^{235}U$ in loop')
|
||||
plt.savefig(f'iso_temperature_feedback_{str(u235_load)}kg.png')
|
||||
|
||||
def rod_bank(temp, positions, loads):
|
||||
load_results = {}
|
||||
for load in loads:
|
||||
y_data = []
|
||||
y_err = []
|
||||
# Rods position in inches where 0 is fully inserted and 51 fully withdrawn
|
||||
for pos in positions:
|
||||
model = build(make_tally=False, plot_geom=False, fuel_temp=temp, u235_load=load, cr1_pos=pos, cr2_pos=pos, cr3_pos=pos)
|
||||
res=model.run()
|
||||
with openmc.StatePoint(res) as sp:
|
||||
y_data.append(sp.keff.n)
|
||||
y_err.append(sp.keff.s)
|
||||
load_results[load] = y_data,y_err
|
||||
df=pd.DataFrame(load_results)
|
||||
df.to_csv('rod_bank')
|
||||
|
||||
def mass_reactivity(loads):
|
||||
y_data = []
|
||||
for load in loads:
|
||||
model = build(make_tally=False, plot_geom=False, fuel_temp=648.9, u235_load=load, cr1_pos=51, cr2_pos=51, cr3_pos=51)
|
||||
res=model.run()
|
||||
with openmc.StatePoint(res) as sp:
|
||||
y_data.append(sp.keff.n)
|
||||
loads = np.array(loads)
|
||||
y_data = np.array(y_data)
|
||||
print(y_data)
|
||||
x_data=(loads-loads[0])/loads[0]
|
||||
y_data=(y_data-y_data[0])/(y_data*y_data[0])
|
||||
popt, pcov = curve_fit(func, x_data, y_data)
|
||||
plt.plot(x_data,y_data,'o',label='data')
|
||||
plt.plot(x_data, func(x_data, *popt), '--', label='fit: f(x)={:.3f}x + {:.3f}'.format(popt[0],popt[1]))
|
||||
plt.legend()
|
||||
plt.ylabel(r'Change of reactivity, $\Delta k /(k_1k_2)$')
|
||||
plt.xlabel(r'Fractional change of $^{235}U$ loading in loop, $\Delta m /m$')
|
||||
plt.savefig('mass_reactivity')
|
||||
|
||||
def triton_adder(mass):
|
||||
|
||||
df = pd.read_csv('ADDER/msre_simple_fuel_mod.csv')
|
||||
|
|
@ -486,7 +609,14 @@ def triton_adder(mass):
|
|||
if __name__ == '__main__':
|
||||
mass = 4590 #tot fuel salt mass [kg]
|
||||
power = 8e6 #total thermal power [W]
|
||||
#run(build(make_tally=True, plot_geom=True), mass, power)
|
||||
control_rod_worth(build(make_tally=False, plot_geom=False))
|
||||
#depletion(build(make_tally=False, plot_geom=False), mass, power)
|
||||
fuel_temp = 638.3#fuel temperature at the initial criticality point reported [C]
|
||||
#run(build(make_tally=False, plot_geom=True, u235_load=71.71, cr1_pos=0,fuel_temp=648.9), mass, power)
|
||||
#control_rod_worth(build(make_tally=False, plot_geom=False, fuel_temp=fuel_temp))
|
||||
#rod_bank(648.9, [0,51],[67.94,69.94,71.71])
|
||||
feedback_isothermal([598.9,648.9,698.9],67.86)
|
||||
feedback_isothermal([598.9,648.9,698.9],71.71)
|
||||
feedback_isothermal([598.9,648.9,698.9],69.85)
|
||||
#feedback_fuel([598.9])#,648.9,698.9])
|
||||
#mass_reactivity([65.25,66,67,68,69,70,71,72])
|
||||
#depletion(build(make_tally=False, plot_geom=False, fuel_temp=fuel_temp), mass, power)
|
||||
#triton_adder(mass)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue