Add new h5m files
This commit is contained in:
church89
parent
1e06779e09
commit
23e0706324
11 changed files with 10 additions and 20890 deletions
Binary file not shown.
BIN
h5m/msre_full.h5m
(Stored with Git LFS)
Normal file
BIN
h5m/msre_full.h5m
(Stored with Git LFS)
Normal file
Binary file not shown.
318
msre_NEW.py
318
msre_NEW.py
|
|
@ -1,318 +0,0 @@
|
|||
import os
|
||||
import openmc
|
||||
import openmc.deplete
|
||||
import openmc.lib
|
||||
import numpy as np
|
||||
from math import log10, sqrt
|
||||
from collections import OrderedDict
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
import math
|
||||
import numpy as np
|
||||
import re
|
||||
from pathlib import Path
|
||||
import h5py
|
||||
|
||||
salt_temp = 648.9
|
||||
salt = openmc.Material(name="salt", temperature = salt_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', 2.32151)
|
||||
#salt.volume = 4560/2.32151*1000
|
||||
#salt.volume = 70 * 0.0283168 *1e6 #Circulating primary salt: 70 ft3 (ORNL-4658)
|
||||
salt.volume = 1.996 * 1e6 # https://info.ornl.gov/sites/publications/Files/Pub173113.pdf
|
||||
#moderator blocks170
|
||||
graphite = openmc.Material(name='graphite',temperature= salt_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= salt_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')
|
||||
inor.add_element('Cr',(6+8)/2,'wo')
|
||||
inor.add_element('Fe',5,'wo')
|
||||
inor.add_element('C',(0.04+0.08)/2,'wo')
|
||||
inor.add_element('Al',0.25,'wo')
|
||||
inor.add_element('Ti',0.25,'wo')
|
||||
inor.add_element('S',0.02,'wo')
|
||||
inor.add_element('Mn',1.0,'wo')
|
||||
inor.add_element('Si',1.0,'wo')
|
||||
inor.add_element('Cu',0.35,'wo')
|
||||
inor.add_element('B',0.010,'wo')
|
||||
inor.add_element('W',0.5,'wo')
|
||||
inor.add_element('P',0.015,'wo')
|
||||
inor.add_element('Co',0.2,'wo')
|
||||
|
||||
#helium
|
||||
helium = openmc.Material(name='helium')
|
||||
helium.add_element('He',1.0)
|
||||
helium.set_density('g/cm3',1.03*(10**-4))
|
||||
|
||||
#Control rods inconel clad
|
||||
trace = 0.01
|
||||
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')
|
||||
inconel.add_element('Mn',0.25,percent_type='wo')
|
||||
inconel.add_element('Si',0.25,percent_type='wo')
|
||||
inconel.add_element('Cu',0.2,percent_type='wo')
|
||||
inconel.add_element('Co',0.2,percent_type='wo')
|
||||
inconel.add_element('Al',0.2,percent_type='wo')
|
||||
inconel.add_element('Ti',0.2,percent_type='wo')
|
||||
inconel.add_element('Ta',0.5,percent_type='wo')
|
||||
inconel.add_element('W',0.5,percent_type='wo')
|
||||
inconel.add_element('Zn',0.2,percent_type='wo')
|
||||
inconel.add_element('Zr',0.1,percent_type='wo')
|
||||
inconel.add_element('C',trace,percent_type='wo')
|
||||
inconel.add_element('Mo',trace,percent_type='wo')
|
||||
inconel.add_element('Ag',trace,percent_type='wo')
|
||||
inconel.add_element('B',trace,percent_type='wo')
|
||||
inconel.add_element('Ba',trace,percent_type='wo')
|
||||
inconel.add_element('Be',trace,percent_type='wo')
|
||||
inconel.add_element('Ca',trace,percent_type='wo')
|
||||
inconel.add_element('Cd',trace,percent_type='wo')
|
||||
inconel.add_element('V',trace,percent_type='wo')
|
||||
inconel.add_element('Sn',trace,percent_type='wo')
|
||||
inconel.add_element('Mg',trace,percent_type='wo')
|
||||
inconel.set_density('g/cm3',8.5)
|
||||
|
||||
# SS 316 control rod flexible hose
|
||||
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')
|
||||
ss316.add_element('Cr',16.55,'wo')
|
||||
ss316.add_element('Cu',0.16,'wo')
|
||||
ss316.add_element('Ni',10,'wo')
|
||||
ss316.add_element('P',0.029,'wo')
|
||||
ss316.add_element('S',0.027,'wo')
|
||||
ss316.add_element('Mo',2.02,'wo')
|
||||
ss316.add_element('N',0.036,'wo')
|
||||
ss316.add_element('Fe',70.622,'wo')
|
||||
ss316.set_density('g/cm3',7.99)
|
||||
|
||||
#Control rods bushing posion material
|
||||
Gd2O3 = openmc.Material()
|
||||
Gd2O3.add_element('Gd',2)
|
||||
Gd2O3.add_element('O',3)
|
||||
Gd2O3.set_density('g/cm3',7.41)
|
||||
Al2O3 = openmc.Material()
|
||||
Al2O3.add_element('Al',2)
|
||||
Al2O3.add_element('O',3)
|
||||
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')
|
||||
concrete.add_element('H',0.005,'wo')
|
||||
concrete.add_element('O',0.496,'wo')
|
||||
concrete.add_element('Si',0.314,'wo')
|
||||
concrete.add_element('Ca',0.083,'wo')
|
||||
concrete.add_element('Na',0.017,'wo')
|
||||
concrete.add_element('Mn',0.002,'wo')
|
||||
concrete.add_element('Al',0.046,'wo')
|
||||
concrete.add_element('S',0.001,'wo')
|
||||
concrete.add_element('K',0.019,'wo')
|
||||
concrete.add_element('Fe',0.012,'wo')
|
||||
concrete.set_density('g/cm3',2.35)
|
||||
|
||||
#Thermal shielding as water and SS305 (50-50)
|
||||
water = openmc.Material()
|
||||
water.add_element('H',2)
|
||||
water.add_element('O',1)
|
||||
water.set_density('g/cm3',0.997)
|
||||
|
||||
#stainless steel 304
|
||||
ss304 = openmc.Material()
|
||||
ss304.add_element('C',0.08,'wo')
|
||||
ss304.add_element('Mn',2,'wo')
|
||||
ss304.add_element('P',0.045,'wo')
|
||||
ss304.add_element('S',0.03,'wo')
|
||||
ss304.add_element('Si',0.75,'wo')
|
||||
ss304.add_element('Cr',19,'wo')
|
||||
ss304.add_element('Ni',10,'wo')
|
||||
ss304.add_element('N',0.1,'wo')
|
||||
ss304.add_element('Fe',67.995, 'wo')
|
||||
ss304.set_density('g/cm3',7.93)
|
||||
shield = openmc.Material.mix_materials([water,ss304],[0.5,0.5],'vo')
|
||||
shield.temperature = 32.2 + 273.15
|
||||
shield.name='steelwater'
|
||||
|
||||
# "Careytemp 1600" by Philip Carey Manufacturing Compamy (Cincinnati)
|
||||
# http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf
|
||||
insulation=openmc.Material(name='insulation')
|
||||
insulation.add_element('Si',1)
|
||||
insulation.add_element('O',2)
|
||||
insulation.set_density('g/cm3',0.16) #https://www.osti.gov/servlets/purl/1411211
|
||||
|
||||
# sand water, not sure about this material
|
||||
sandwater=openmc.Material(name='sandwater')
|
||||
sandwater.add_element('Fe',3)
|
||||
sandwater.add_element('O',4)
|
||||
sandwater.set_density('g/cm3',6)
|
||||
|
||||
#Vessel anular steel
|
||||
steel = openmc.Material(name='steel')
|
||||
steel.add_element('Fe',1)
|
||||
steel.set_density('g/cm3',7.85)
|
||||
|
||||
mats = openmc.Materials([salt, graphite, inor, helium, inconel, shield, concrete,
|
||||
steel, ss316, sandwater, insulation, bush])
|
||||
|
||||
# CAD h5m files
|
||||
core_h5m = 'h5m/msre_reactor_1e-2.h5m'
|
||||
control_rod1_h5m = 'h5m/msre_control_rod_1e-2.h5m'
|
||||
|
||||
#Geometry
|
||||
core = openmc.DAGMCUniverse(filename = core_h5m, auto_geom_ids = True,
|
||||
universe_id = 1)
|
||||
control_rod1 = openmc.DAGMCUniverse(filename = control_rod1_h5m,
|
||||
auto_geom_ids = True, universe_id=2)
|
||||
core_region = core.bounding_region()
|
||||
cr1_region = control_rod1.bounding_region(boundary_type = 'transmission',
|
||||
starting_id=20000)
|
||||
|
||||
# Extend control rod region1 to include upwards translations
|
||||
cr1_region = cr1_region | cr1_region.translate([0, 0, 150])
|
||||
|
||||
# Shift control rods 1 by offset to defined control rod 2 and 3 regions
|
||||
offset = 10.163255 #cm, offset between cr1, cr2 and cr3
|
||||
cr2_region = cr1_region.translate([-offset, 0, 0])
|
||||
cr3_region = cr1_region.translate([-offset, offset, 0])
|
||||
|
||||
#Define cells
|
||||
core_cell = openmc.Cell(region=~(cr1_region | cr2_region | cr3_region) & core_region ,
|
||||
fill=core)
|
||||
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)
|
||||
|
||||
#Fix control rods initial positions
|
||||
start_pos = 19.2 #cm, geometrical distance between lower bottom and starting point
|
||||
top_pos = 51 * 2.54 # cm, initial position of control rod1 with respect to start post
|
||||
init_pos = 40 * 2.54
|
||||
setattr(cr1_cell, 'translation', [0, 0, start_pos + init_pos])
|
||||
setattr(cr2_cell, 'translation', [-offset, 0, start_pos + top_pos])
|
||||
setattr(cr3_cell, 'translation', [-offset, offset, start_pos + top_pos])
|
||||
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 = 80
|
||||
settings.inactive = 20
|
||||
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)
|
||||
|
||||
#Tallies
|
||||
tally = openmc.Tally(name="heating")
|
||||
if settings.photon_transport:
|
||||
heating_score = 'heating'
|
||||
else:
|
||||
heating_score = 'heating-local'
|
||||
tally.scores.append(heating_score)
|
||||
tallies = openmc.Tallies([tally])
|
||||
|
||||
# Depletion settings
|
||||
model = openmc.model.Model(geometry,mats,settings,tallies)
|
||||
|
||||
#Plots
|
||||
colors = {salt:'yellow', graphite:'black', inor: 'grey', helium: 'cyan', inconel: 'grey',
|
||||
bush: 'blue', ss316: 'grey', concrete: 'brown', shield: 'red', insulation: 'green',
|
||||
sandwater: 'lightgreen', steel: 'grey'}
|
||||
plots = openmc.Plots()
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xy'
|
||||
plot.width = (150,150)
|
||||
plot.pixels = (1500,1500)
|
||||
plot.origin = (0,0,160)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (750,1500)
|
||||
plot.origin = (0,-5,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'yz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (500,1000)
|
||||
plot.origin = (5,0,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
|
||||
op = openmc.deplete.CoupledOperator(model,
|
||||
normalization_mode = "energy-deposition")
|
||||
|
||||
timesteps = [5,5,30,30,30,180,95]
|
||||
power = 8 #7.34
|
||||
integrator = openmc.deplete.CECMIntegrator(op, timesteps=timesteps,
|
||||
timestep_units='d', power=power)
|
||||
|
||||
integrator.add_transfer_rate('salt', ['Xe','Kr'], 4.067e-5)
|
||||
integrator.add_transfer_rate('salt', ['Se','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Sb','Te'], 8.777e-3)
|
||||
integrator.add_redox('salt', {'Be9':1})
|
||||
# integrator.add_batchwise( 'CR1', 'translation', axis = 2,
|
||||
# bracket = [-2, 5], #cm
|
||||
# bracket_limit = [0, top_pos+start_pos],
|
||||
# #atom_density_limit = 1e-14, #atoms/b-cm
|
||||
# tol = 0.1)
|
||||
|
||||
|
||||
# integrator.add_batchwise('refuel', mats_id_or_name = ['salt'],
|
||||
# mat_vector = {'U235': 1},
|
||||
# bracket = [1e2,1e3], #grams
|
||||
# bracket_limit = [0,1e5], #grams
|
||||
# tol = 0.01,
|
||||
# restart_level = start_pos + init_pos)
|
||||
|
||||
#integrator.add_batchwise_wrap('1', interrupt=True)
|
||||
|
||||
integrator.integrate()
|
||||
BIN
msre_control_rod.h5m
(Stored with Git LFS)
Normal file
BIN
msre_control_rod.h5m
(Stored with Git LFS)
Normal file
Binary file not shown.
|
|
@ -1,194 +0,0 @@
|
|||
import openmc
|
||||
import openmc.deplete
|
||||
import os
|
||||
import re
|
||||
import matplotlib.pyplot as plt
|
||||
import seaborn as sns
|
||||
import datetime
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
from pathlib import Path
|
||||
regex = re.compile(r'(\d+|\s+)')
|
||||
|
||||
path_to_results = '/home/lorenzo/mnt/uranium/msre/depletion_results.h5'
|
||||
save_dir = Path(os.path.realpath(path_to_results)).parent
|
||||
mats = {mat.name: mat.id for mat in openmc.material.Materials.from_xml(save_dir / 'materials.xml')
|
||||
if mat.depletable}
|
||||
|
||||
results = openmc.deplete.Results(path_to_results)
|
||||
_, keff = results.get_keff()
|
||||
n_xe = 0
|
||||
n_kr = 0
|
||||
for nuc,_ in openmc.data.isotopes('Xe'):
|
||||
n_xe += results.get_atoms(str(mats['salt']), nuc)[1]
|
||||
for nuc,_ in openmc.data.isotopes('Kr'):
|
||||
n_kr += results.get_atoms(str(mats['salt']), nuc)[1]
|
||||
|
||||
df=pd.read_excel('/home/lorenzo/mnt/uranium/msre/MSRE_235_233_Power_History_R24E.xlsx')
|
||||
# U235 operation run
|
||||
dt = df['Duration (h)'][:94]
|
||||
# First operation in MW range(ORNL-4674)
|
||||
t0=datetime.datetime.strptime('24/01/1966', '%d/%m/%Y')
|
||||
date_times = []
|
||||
for t in dt:
|
||||
date_times.append(t0)
|
||||
t0 += datetime.timedelta(hours=t)
|
||||
date_times.append(t0)
|
||||
|
||||
power=df['Power (MWth)'][:95]
|
||||
plt.figure(figsize=(18,10))
|
||||
ax = plt.subplot()
|
||||
ax1 = ax.twinx()
|
||||
ax.errorbar(date_times, [k[0] for k in keff], [k[1] for k in keff], marker='o', color='black', markersize=4, fmt=' ')
|
||||
|
||||
ax1.step(date_times, power, where='post', color='red')
|
||||
ax.set_ylabel(r'$k_{eff}\pm\sigma$')
|
||||
ax.set_ylim(0.98,1.01)
|
||||
ax.tick_params(axis='y')
|
||||
ax1.set_ylabel('Power [MWth]',color='red')
|
||||
ax1.tick_params(axis='y', colors='red')
|
||||
plt.savefig(f'{save_dir}/keff', dpi=600)
|
||||
|
||||
# Microscopic absorption cross section at 0.0253 eV
|
||||
xs_xe135 = 2664214.0
|
||||
xs_u235 = 686.006994850397
|
||||
_, n_xe135 = results.get_atoms(str(mats['salt']), 'Xe135')
|
||||
_, n_u235 = results.get_atoms(str(mats['salt']), 'U235')
|
||||
# Poison fraction
|
||||
pf = (xs_xe135*n_xe135)/(xs_u235*n_u235)*100
|
||||
plt.figure(figsize=(18,10))
|
||||
plt.plot(date_times, pf)
|
||||
plt.ylabel('Xe posion fraction [%]')
|
||||
plt.savefig(f'{save_dir}/fission_products', dpi=600)
|
||||
|
||||
inventory = dict()
|
||||
for nuc,_ in openmc.data.isotopes('U'):
|
||||
inventory[nuc] = results.get_atoms(str(mats['salt']), nuc)[1] / openmc.data.AVOGADRO * openmc.data.atomic_mass(nuc) / 1000
|
||||
|
||||
for nuc in ['Pu238','Pu239','Pu240','Pu241','Pu242']:
|
||||
inventory[nuc] = results.get_atoms(str(mats['salt']), nuc)[1] / openmc.data.AVOGADRO * openmc.data.atomic_mass(nuc) / 1000
|
||||
|
||||
plt.figure(figsize=(18,10))
|
||||
for nuc, mass in inventory.items():
|
||||
plt.plot(date_times, mass, label=nuc)
|
||||
plt.ylabel('Mass [g]')
|
||||
plt.yscale('log')
|
||||
plt.ylim(1e-5)
|
||||
plt.legend()
|
||||
plt.savefig(f'{save_dir}/inventory', dpi=600)
|
||||
|
||||
lib = openmc.data.DataLibrary.from_xml(os.environ.get("OPENMC_CROSS_SECTION"))
|
||||
nuclides = set()
|
||||
for library in lib.libraries:
|
||||
if library['type'] != 'neutron':
|
||||
continue
|
||||
for name in library['materials']:
|
||||
if name not in nuclides:
|
||||
nuclides.add(name)
|
||||
|
||||
# Let's begin by making some useful groupings
|
||||
gaseos = ['H', 'He', 'Ne', 'Ar', 'Kr', 'Xe', 'Rn'] #gaseous fission products
|
||||
noble_metals = ['Se','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Sb','Te'] # noble metals fission products
|
||||
metals = ['Cr','Mn','Fe','Co','Ni','Cu','Zn','Hf','Zr','W',]
|
||||
halogens = ['F','Cl','Br','I','At']
|
||||
alkali_metals = ['Li','Na','K','Rb','Cs']
|
||||
alkali_earths= ['Be','Mg','Ca','Sr','Ba','Ra']
|
||||
lanthanides = ['Y','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu']
|
||||
m_a = ['Ac','Th','Pa','Np','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr']
|
||||
|
||||
# Get fissile nuclides in the fuel, based on Ronen's rule for determining fissile isotopes
|
||||
fissile = []
|
||||
|
||||
for nuc in nuclides:
|
||||
elm = regex.split(nuc)[0]
|
||||
a = round(openmc.data.atomic_mass(nuc))
|
||||
z = openmc.data.ATOMIC_NUMBER[elm]
|
||||
if 90 <= z <= 100:
|
||||
ronen = 2*z -(a-z)
|
||||
if ronen in [41,43,45]:
|
||||
fissile.append(nuc)
|
||||
|
||||
# Calculate totat absorption rate of fissile nuclides
|
||||
tot_abs_rate = 0
|
||||
for nuc in fissile:
|
||||
tot_abs_rate += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
tot_abs_rate += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
|
||||
nuclides_stack = dict()
|
||||
groups_stack = {'Gaseos':0, 'Noble metals':0, 'Metals':0, 'Halogens':0 , 'Alkali metals':0, 'Alkali earths':0, 'Lanthanides':0, 'MA':0, 'Others':0}
|
||||
for nuc in nuclides:
|
||||
|
||||
if regex.split(nuc)[0] in ['U','Pu']:
|
||||
nuclides_stack[nuc] = results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
nuclides_stack[nuc] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
nuclides_stack[nuc] /= tot_abs_rate
|
||||
|
||||
|
||||
elif regex.split(nuc)[0] in gaseos:
|
||||
groups_stack['Gaseos'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Gaseos'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in noble_metals:
|
||||
groups_stack['Noble metals'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Noble metals'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in metals:
|
||||
groups_stack['Metals'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Metals'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in halogens:
|
||||
groups_stack['Halogens'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Halogens'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in alkali_metals:
|
||||
groups_stack['Alkali metals'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Alkali metals'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in alkali_earths:
|
||||
groups_stack['Alkali earths'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Alkali earths'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in lanthanides:
|
||||
groups_stack['Lanthanides'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Lanthanides'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
elif regex.split(nuc)[0] in m_a:
|
||||
groups_stack['MA'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['MA'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
else:
|
||||
groups_stack['Others'] += results.get_reaction_rate(str(mats['salt']), nuc, 'fission')[1]
|
||||
groups_stack['Others'] += results.get_reaction_rate(str(mats['salt']), nuc, '(n,gamma)')[1]
|
||||
|
||||
# Divide each array by the total absorption reaction rate of fissile nuclides
|
||||
for g in groups_stack.keys():
|
||||
groups_stack[g] /= tot_abs_rate
|
||||
|
||||
# Sort dictionary groups
|
||||
groups_stack=dict(reversed(sorted(groups_stack.items(), key=lambda item: item[1][len(item)])))
|
||||
|
||||
# Create red color palette for groups_stack
|
||||
colors = list(reversed(sns.color_palette("Reds", len(groups_stack))))
|
||||
|
||||
# Order uranium series
|
||||
u_series = {key:value for key,value in nuclides_stack.items() if key.startswith('U')}
|
||||
u_series = dict(reversed(sorted(u_series.items(), key=lambda item: item[1][len(item)])))
|
||||
# Create green color palette for Uranium isotopes
|
||||
colors += list(reversed(sns.color_palette("Greens", len(u_series))))
|
||||
|
||||
# Order plutionium series
|
||||
pu_series = {key:value for key,value in nuclides_stack.items() if key.startswith('Pu')}
|
||||
pu_series = dict(reversed(sorted(pu_series.items(), key=lambda item: item[1][len(item)])))
|
||||
# Create blue color palette for plutonium isotopes
|
||||
colors += list(reversed(sns.color_palette("Blues", len(pu_series))))
|
||||
# Add uramium and plutonium series to the stack
|
||||
groups_stack.update(u_series)
|
||||
groups_stack.update(pu_series)
|
||||
|
||||
plt.figure(figsize=(15,10))
|
||||
plt.stackplot(date_times, groups_stack.values(), labels=groups_stack.keys(),
|
||||
edgecolor="black", linewidth=0.5,colors=colors, alpha=0.8)
|
||||
handles, labels = plt.gca().get_legend_handles_labels()
|
||||
legend = plt.legend([handles[idx] for idx in list(reversed(np.arange(0,len(handles),1)))],
|
||||
[labels[idx] for idx in list(reversed(np.arange(0,len(handles),1)))],
|
||||
bbox_to_anchor=(1.05,1), loc='upper left',
|
||||
borderaxespad=0, ncol=2,
|
||||
fontsize=15)
|
||||
plt.title('Fuel salt neutrons absorption distribution per neutron absorbed in fissile isotopes',
|
||||
weight='bold', fontsize=17)
|
||||
plt.xticks(fontsize=13)
|
||||
plt.yticks(fontsize=13)
|
||||
plt.tight_layout()
|
||||
plt.savefig(f'{save_dir}/neutrons', dpi=600)
|
||||
|
|
@ -1,339 +0,0 @@
|
|||
import os
|
||||
import openmc
|
||||
import openmc.deplete
|
||||
import openmc.lib
|
||||
import numpy as np
|
||||
from math import log10, sqrt
|
||||
from collections import OrderedDict
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
import math
|
||||
import numpy as np
|
||||
import re
|
||||
|
||||
salt_temp = 648.9
|
||||
salt = openmc.Material(name="salt", temperature = salt_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', 2.32151)
|
||||
#salt.volume = 4560/2.32151*1000
|
||||
#salt.volume = 70 * 0.0283168 *1e6 #Circulating primary salt: 70 ft3 (ORNL-4658)
|
||||
salt.volume = 1.996 * 1e6 # https://info.ornl.gov/sites/publications/Files/Pub173113.pdf
|
||||
#moderator blocks170
|
||||
graphite = openmc.Material(name='graphite',temperature= salt_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= salt_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')
|
||||
inor.add_element('Cr',(6+8)/2,'wo')
|
||||
inor.add_element('Fe',5,'wo')
|
||||
inor.add_element('C',(0.04+0.08)/2,'wo')
|
||||
inor.add_element('Al',0.25,'wo')
|
||||
inor.add_element('Ti',0.25,'wo')
|
||||
inor.add_element('S',0.02,'wo')
|
||||
inor.add_element('Mn',1.0,'wo')
|
||||
inor.add_element('Si',1.0,'wo')
|
||||
inor.add_element('Cu',0.35,'wo')
|
||||
inor.add_element('B',0.010,'wo')
|
||||
inor.add_element('W',0.5,'wo')
|
||||
inor.add_element('P',0.015,'wo')
|
||||
inor.add_element('Co',0.2,'wo')
|
||||
|
||||
#helium
|
||||
helium = openmc.Material(name='helium')
|
||||
helium.add_element('He',1.0)
|
||||
helium.set_density('g/cm3',1.03*(10**-4))
|
||||
|
||||
#Control rods inconel clad
|
||||
trace = 0.01
|
||||
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')
|
||||
inconel.add_element('Mn',0.25,percent_type='wo')
|
||||
inconel.add_element('Si',0.25,percent_type='wo')
|
||||
inconel.add_element('Cu',0.2,percent_type='wo')
|
||||
inconel.add_element('Co',0.2,percent_type='wo')
|
||||
inconel.add_element('Al',0.2,percent_type='wo')
|
||||
inconel.add_element('Ti',0.2,percent_type='wo')
|
||||
inconel.add_element('Ta',0.5,percent_type='wo')
|
||||
inconel.add_element('W',0.5,percent_type='wo')
|
||||
inconel.add_element('Zn',0.2,percent_type='wo')
|
||||
inconel.add_element('Zr',0.1,percent_type='wo')
|
||||
inconel.add_element('C',trace,percent_type='wo')
|
||||
inconel.add_element('Mo',trace,percent_type='wo')
|
||||
inconel.add_element('Ag',trace,percent_type='wo')
|
||||
inconel.add_element('B',trace,percent_type='wo')
|
||||
inconel.add_element('Ba',trace,percent_type='wo')
|
||||
inconel.add_element('Be',trace,percent_type='wo')
|
||||
inconel.add_element('Ca',trace,percent_type='wo')
|
||||
inconel.add_element('Cd',trace,percent_type='wo')
|
||||
inconel.add_element('V',trace,percent_type='wo')
|
||||
inconel.add_element('Sn',trace,percent_type='wo')
|
||||
inconel.add_element('Mg',trace,percent_type='wo')
|
||||
inconel.set_density('g/cm3',8.5)
|
||||
|
||||
# SS 316 control rod flexible hose
|
||||
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')
|
||||
ss316.add_element('Cr',16.55,'wo')
|
||||
ss316.add_element('Cu',0.16,'wo')
|
||||
ss316.add_element('Ni',10,'wo')
|
||||
ss316.add_element('P',0.029,'wo')
|
||||
ss316.add_element('S',0.027,'wo')
|
||||
ss316.add_element('Mo',2.02,'wo')
|
||||
ss316.add_element('N',0.036,'wo')
|
||||
ss316.add_element('Fe',70.622,'wo')
|
||||
ss316.set_density('g/cm3',7.99)
|
||||
|
||||
#Control rods bushing posion material
|
||||
Gd2O3 = openmc.Material()
|
||||
Gd2O3.add_element('Gd',2)
|
||||
Gd2O3.add_element('O',3)
|
||||
Gd2O3.set_density('g/cm3',7.41)
|
||||
Al2O3 = openmc.Material()
|
||||
Al2O3.add_element('Al',2)
|
||||
Al2O3.add_element('O',3)
|
||||
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')
|
||||
concrete.add_element('H',0.005,'wo')
|
||||
concrete.add_element('O',0.496,'wo')
|
||||
concrete.add_element('Si',0.314,'wo')
|
||||
concrete.add_element('Ca',0.083,'wo')
|
||||
concrete.add_element('Na',0.017,'wo')
|
||||
concrete.add_element('Mn',0.002,'wo')
|
||||
concrete.add_element('Al',0.046,'wo')
|
||||
concrete.add_element('S',0.001,'wo')
|
||||
concrete.add_element('K',0.019,'wo')
|
||||
concrete.add_element('Fe',0.012,'wo')
|
||||
concrete.set_density('g/cm3',2.35)
|
||||
|
||||
#Thermal shielding as water and SS305 (50-50)
|
||||
water = openmc.Material()
|
||||
water.add_element('H',2)
|
||||
water.add_element('O',1)
|
||||
water.set_density('g/cm3',0.997)
|
||||
|
||||
#stainless steel 304
|
||||
ss304 = openmc.Material()
|
||||
ss304.add_element('C',0.08,'wo')
|
||||
ss304.add_element('Mn',2,'wo')
|
||||
ss304.add_element('P',0.045,'wo')
|
||||
ss304.add_element('S',0.03,'wo')
|
||||
ss304.add_element('Si',0.75,'wo')
|
||||
ss304.add_element('Cr',19,'wo')
|
||||
ss304.add_element('Ni',10,'wo')
|
||||
ss304.add_element('N',0.1,'wo')
|
||||
ss304.add_element('Fe',67.995, 'wo')
|
||||
ss304.set_density('g/cm3',7.93)
|
||||
shield = openmc.Material.mix_materials([water,ss304],[0.5,0.5],'vo')
|
||||
shield.temperature = 32.2 + 273.15
|
||||
shield.name='steelwater'
|
||||
|
||||
# "Careytemp 1600" by Philip Carey Manufacturing Compamy (Cincinnati)
|
||||
# http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf
|
||||
insulation=openmc.Material(name='insulation')
|
||||
insulation.add_element('Si',1)
|
||||
insulation.add_element('O',2)
|
||||
insulation.set_density('g/cm3',0.16) #https://www.osti.gov/servlets/purl/1411211
|
||||
|
||||
# sand water, not sure about this material
|
||||
sandwater=openmc.Material(name='sandwater')
|
||||
sandwater.add_element('Fe',3)
|
||||
sandwater.add_element('O',4)
|
||||
sandwater.set_density('g/cm3',6)
|
||||
|
||||
#Vessel anular steel
|
||||
steel = openmc.Material(name='steel')
|
||||
steel.add_element('Fe',1)
|
||||
steel.set_density('g/cm3',7.85)
|
||||
|
||||
mats = openmc.Materials([salt, graphite, inor, helium, inconel, shield, concrete,
|
||||
steel, ss316, sandwater, insulation, bush])
|
||||
|
||||
|
||||
# CAD h5m files
|
||||
core_h5m = 'h5m/msre_reactor_1e-2.h5m'
|
||||
control_rod1_h5m = 'h5m/msre_control_rod_1e-2.h5m'
|
||||
|
||||
#Geometry
|
||||
core = openmc.DAGMCUniverse(filename = core_h5m, auto_geom_ids = True,
|
||||
universe_id = 1)
|
||||
control_rod1 = openmc.DAGMCUniverse(filename = control_rod1_h5m,
|
||||
auto_geom_ids = True, universe_id=2)
|
||||
core_region = core.bounding_region()
|
||||
cr1_region = control_rod1.bounding_region(boundary_type = 'transmission',
|
||||
starting_id=20000)
|
||||
|
||||
# Extend control rod region1 to include upwards translations
|
||||
cr1_region = cr1_region | cr1_region.translate([0, 0, 150])
|
||||
|
||||
# Shift control rods 1 by offset to defined control rod 2 and 3 regions
|
||||
offset = 10.163255 #cm, offset between cr1, cr2 and cr3
|
||||
cr2_region = cr1_region.translate([-offset, 0, 0])
|
||||
cr3_region = cr1_region.translate([-offset, offset, 0])
|
||||
|
||||
#Define cells
|
||||
core_cell = openmc.Cell(region=~(cr1_region | cr2_region | cr3_region) & core_region ,
|
||||
fill=core)
|
||||
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)
|
||||
|
||||
#Fix control rods initial positions
|
||||
start_pos = 19.2 #cm, geometrical distance between lower bottom and starting point
|
||||
top_pos = 51 * 2.54 # cm, initial position of control rod1 with respect to start post
|
||||
init_pos = 35 * 2.54
|
||||
setattr(cr1_cell, 'translation', [0, 0, start_pos + init_pos])
|
||||
setattr(cr2_cell, 'translation', [-offset, 0, start_pos + top_pos])
|
||||
setattr(cr3_cell, 'translation', [-offset, offset, start_pos + top_pos])
|
||||
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 = 60
|
||||
settings.inactive = 20
|
||||
settings.particles = 500000
|
||||
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)
|
||||
|
||||
#Tallies
|
||||
tally = openmc.Tally(name="heating")
|
||||
if settings.photon_transport:
|
||||
heating_score = 'heating'
|
||||
else:
|
||||
heating_score = 'heating-local'
|
||||
tally.scores.append(heating_score)
|
||||
tallies = openmc.Tallies([tally])
|
||||
|
||||
# Depletion settings
|
||||
model = openmc.model.Model(geometry,mats,settings,tallies)
|
||||
#model.export_to_xml()
|
||||
|
||||
#Plots
|
||||
colors = {salt:'yellow', graphite:'black', inor: 'grey', helium: 'cyan', inconel: 'grey',
|
||||
bush: 'blue', ss316: 'grey', concrete: 'brown', shield: 'red', insulation: 'green',
|
||||
sandwater: 'lightgreen', steel: 'grey'}
|
||||
plots = openmc.Plots()
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xy'
|
||||
plot.width = (150,150)
|
||||
plot.pixels = (1500,1500)
|
||||
plot.origin = (0,0,160)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (750,1500)
|
||||
plot.origin = (0,-5,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'yz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (500,1000)
|
||||
plot.origin = (5,0,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
#model.plot_geometry()
|
||||
|
||||
#results=model.run()
|
||||
|
||||
op = openmc.deplete.CoupledOperator(model,
|
||||
normalization_mode = "energy-deposition")
|
||||
|
||||
df=pd.read_excel('MSRE_235_233_Power_History_R24E.xlsx')
|
||||
|
||||
#U235
|
||||
timesteps = df['Duration (d)'][:94].values
|
||||
power = df['Power (MWth)'][:94].values*1000000
|
||||
|
||||
#Add one further timestep at the end of every power run
|
||||
# timesteps_ext = []
|
||||
# power_ext = []
|
||||
# for i in range(len(timesteps)):
|
||||
# power_ext.append(power[i])
|
||||
# if power[i] != 0.0:
|
||||
# # duplicate power value
|
||||
# power_ext.append(power[i])
|
||||
# # add timestep - 1 sec
|
||||
# timesteps_ext.append(timesteps[i] - 1/3600/24)
|
||||
# # add 1 sec
|
||||
# timesteps_ext.append(1/3600/24)
|
||||
# else:
|
||||
# timesteps_ext.append(timesteps[i])
|
||||
|
||||
integrator = openmc.deplete.CECMIntegrator(op, timesteps=timesteps,
|
||||
timestep_units='d', power=power)
|
||||
|
||||
integrator.add_transfer_rate('salt', ['Xe','Kr'], 4.067e-5)
|
||||
integrator.add_transfer_rate('salt', ['Se','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Sb','Te'], 8.777e-3)
|
||||
|
||||
integrator.add_batchwise('trans', axis = 2, cell_id_or_name = 'CR1',
|
||||
bracket = [-2, 5], #cm
|
||||
bracket_limit = [0, top_pos+start_pos],
|
||||
atom_density_limit = 1e-14, #atoms/b-cm
|
||||
tol = 0.1)
|
||||
|
||||
integrator.add_batchwise('refuel', mats_id_or_name = ['salt'],
|
||||
mat_vector = {'U235': 1},
|
||||
bracket = [1e2,1e3], #grams
|
||||
bracket_limit = [0,1e5], #grams
|
||||
tol = 0.01,
|
||||
restart_level = start_pos + init_pos)
|
||||
|
||||
integrator.add_batchwise_wrap('1')
|
||||
|
||||
integrator.integrate()
|
||||
|
|
@ -1,365 +0,0 @@
|
|||
import os
|
||||
import openmc
|
||||
import openmc.deplete
|
||||
import openmc.lib
|
||||
import numpy as np
|
||||
from math import log10, sqrt
|
||||
from collections import OrderedDict
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
import math
|
||||
import numpy as np
|
||||
import re
|
||||
from pathlib import Path
|
||||
import h5py
|
||||
|
||||
def get_geom_level_from_res(depletion_path, index):
|
||||
|
||||
if depletion_path is None:
|
||||
depletion_path = Path(os.getcwd())
|
||||
|
||||
with h5py.File(depletion_path / 'msr_results.h5','r') as f:
|
||||
items = {}
|
||||
for key in f.keys():
|
||||
if re.split(r'_', key)[0] == 'geometry':
|
||||
items[int(re.split(r'_', key)[1])] = np.array(f.get(key))
|
||||
items = OrderedDict(sorted(items.items()))
|
||||
|
||||
res = list(items.items())[index][1].mean()
|
||||
return res
|
||||
|
||||
salt_temp = 648.9
|
||||
salt = openmc.Material(name="salt", temperature = salt_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', 2.32151)
|
||||
#salt.volume = 4560/2.32151*1000
|
||||
#salt.volume = 70 * 0.0283168 *1e6 #Circulating primary salt: 70 ft3 (ORNL-4658)
|
||||
salt.volume = 1.996 * 1e6 # https://info.ornl.gov/sites/publications/Files/Pub173113.pdf
|
||||
#moderator blocks170
|
||||
graphite = openmc.Material(name='graphite',temperature= salt_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= salt_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')
|
||||
inor.add_element('Cr',(6+8)/2,'wo')
|
||||
inor.add_element('Fe',5,'wo')
|
||||
inor.add_element('C',(0.04+0.08)/2,'wo')
|
||||
inor.add_element('Al',0.25,'wo')
|
||||
inor.add_element('Ti',0.25,'wo')
|
||||
inor.add_element('S',0.02,'wo')
|
||||
inor.add_element('Mn',1.0,'wo')
|
||||
inor.add_element('Si',1.0,'wo')
|
||||
inor.add_element('Cu',0.35,'wo')
|
||||
inor.add_element('B',0.010,'wo')
|
||||
inor.add_element('W',0.5,'wo')
|
||||
inor.add_element('P',0.015,'wo')
|
||||
inor.add_element('Co',0.2,'wo')
|
||||
|
||||
#helium
|
||||
helium = openmc.Material(name='helium')
|
||||
helium.add_element('He',1.0)
|
||||
helium.set_density('g/cm3',1.03*(10**-4))
|
||||
|
||||
#Control rods inconel clad
|
||||
trace = 0.01
|
||||
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')
|
||||
inconel.add_element('Mn',0.25,percent_type='wo')
|
||||
inconel.add_element('Si',0.25,percent_type='wo')
|
||||
inconel.add_element('Cu',0.2,percent_type='wo')
|
||||
inconel.add_element('Co',0.2,percent_type='wo')
|
||||
inconel.add_element('Al',0.2,percent_type='wo')
|
||||
inconel.add_element('Ti',0.2,percent_type='wo')
|
||||
inconel.add_element('Ta',0.5,percent_type='wo')
|
||||
inconel.add_element('W',0.5,percent_type='wo')
|
||||
inconel.add_element('Zn',0.2,percent_type='wo')
|
||||
inconel.add_element('Zr',0.1,percent_type='wo')
|
||||
inconel.add_element('C',trace,percent_type='wo')
|
||||
inconel.add_element('Mo',trace,percent_type='wo')
|
||||
inconel.add_element('Ag',trace,percent_type='wo')
|
||||
inconel.add_element('B',trace,percent_type='wo')
|
||||
inconel.add_element('Ba',trace,percent_type='wo')
|
||||
inconel.add_element('Be',trace,percent_type='wo')
|
||||
inconel.add_element('Ca',trace,percent_type='wo')
|
||||
inconel.add_element('Cd',trace,percent_type='wo')
|
||||
inconel.add_element('V',trace,percent_type='wo')
|
||||
inconel.add_element('Sn',trace,percent_type='wo')
|
||||
inconel.add_element('Mg',trace,percent_type='wo')
|
||||
inconel.set_density('g/cm3',8.5)
|
||||
|
||||
# SS 316 control rod flexible hose
|
||||
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')
|
||||
ss316.add_element('Cr',16.55,'wo')
|
||||
ss316.add_element('Cu',0.16,'wo')
|
||||
ss316.add_element('Ni',10,'wo')
|
||||
ss316.add_element('P',0.029,'wo')
|
||||
ss316.add_element('S',0.027,'wo')
|
||||
ss316.add_element('Mo',2.02,'wo')
|
||||
ss316.add_element('N',0.036,'wo')
|
||||
ss316.add_element('Fe',70.622,'wo')
|
||||
ss316.set_density('g/cm3',7.99)
|
||||
|
||||
#Control rods bushing posion material
|
||||
Gd2O3 = openmc.Material()
|
||||
Gd2O3.add_element('Gd',2)
|
||||
Gd2O3.add_element('O',3)
|
||||
Gd2O3.set_density('g/cm3',7.41)
|
||||
Al2O3 = openmc.Material()
|
||||
Al2O3.add_element('Al',2)
|
||||
Al2O3.add_element('O',3)
|
||||
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')
|
||||
concrete.add_element('H',0.005,'wo')
|
||||
concrete.add_element('O',0.496,'wo')
|
||||
concrete.add_element('Si',0.314,'wo')
|
||||
concrete.add_element('Ca',0.083,'wo')
|
||||
concrete.add_element('Na',0.017,'wo')
|
||||
concrete.add_element('Mn',0.002,'wo')
|
||||
concrete.add_element('Al',0.046,'wo')
|
||||
concrete.add_element('S',0.001,'wo')
|
||||
concrete.add_element('K',0.019,'wo')
|
||||
concrete.add_element('Fe',0.012,'wo')
|
||||
concrete.set_density('g/cm3',2.35)
|
||||
|
||||
#Thermal shielding as water and SS305 (50-50)
|
||||
water = openmc.Material()
|
||||
water.add_element('H',2)
|
||||
water.add_element('O',1)
|
||||
water.set_density('g/cm3',0.997)
|
||||
|
||||
#stainless steel 304
|
||||
ss304 = openmc.Material()
|
||||
ss304.add_element('C',0.08,'wo')
|
||||
ss304.add_element('Mn',2,'wo')
|
||||
ss304.add_element('P',0.045,'wo')
|
||||
ss304.add_element('S',0.03,'wo')
|
||||
ss304.add_element('Si',0.75,'wo')
|
||||
ss304.add_element('Cr',19,'wo')
|
||||
ss304.add_element('Ni',10,'wo')
|
||||
ss304.add_element('N',0.1,'wo')
|
||||
ss304.add_element('Fe',67.995, 'wo')
|
||||
ss304.set_density('g/cm3',7.93)
|
||||
shield = openmc.Material.mix_materials([water,ss304],[0.5,0.5],'vo')
|
||||
shield.temperature = 32.2 + 273.15
|
||||
shield.name='steelwater'
|
||||
|
||||
# "Careytemp 1600" by Philip Carey Manufacturing Compamy (Cincinnati)
|
||||
# http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf
|
||||
insulation=openmc.Material(name='insulation')
|
||||
insulation.add_element('Si',1)
|
||||
insulation.add_element('O',2)
|
||||
insulation.set_density('g/cm3',0.16) #https://www.osti.gov/servlets/purl/1411211
|
||||
|
||||
# sand water, not sure about this material
|
||||
sandwater=openmc.Material(name='sandwater')
|
||||
sandwater.add_element('Fe',3)
|
||||
sandwater.add_element('O',4)
|
||||
sandwater.set_density('g/cm3',6)
|
||||
|
||||
#Vessel anular steel
|
||||
steel = openmc.Material(name='steel')
|
||||
steel.add_element('Fe',1)
|
||||
steel.set_density('g/cm3',7.85)
|
||||
|
||||
mats = openmc.Materials([salt, graphite, inor, helium, inconel, shield, concrete,
|
||||
steel, ss316, sandwater, insulation, bush])
|
||||
|
||||
|
||||
# CAD h5m files
|
||||
core_h5m = 'h5m/msre_reactor_1e-2.h5m'
|
||||
control_rod1_h5m = 'h5m/msre_control_rod_1e-2.h5m'
|
||||
|
||||
#Geometry
|
||||
core = openmc.DAGMCUniverse(filename = core_h5m, auto_geom_ids = True,
|
||||
universe_id = 1)
|
||||
control_rod1 = openmc.DAGMCUniverse(filename = control_rod1_h5m,
|
||||
auto_geom_ids = True, universe_id=2)
|
||||
core_region = core.bounding_region()
|
||||
cr1_region = control_rod1.bounding_region(boundary_type = 'transmission',
|
||||
starting_id=20000)
|
||||
|
||||
# Extend control rod region1 to include upwards translations
|
||||
cr1_region = cr1_region | cr1_region.translate([0, 0, 150])
|
||||
|
||||
# Shift control rods 1 by offset to defined control rod 2 and 3 regions
|
||||
offset = 10.163255 #cm, offset between cr1, cr2 and cr3
|
||||
cr2_region = cr1_region.translate([-offset, 0, 0])
|
||||
cr3_region = cr1_region.translate([-offset, offset, 0])
|
||||
|
||||
#Define cells
|
||||
core_cell = openmc.Cell(region=~(cr1_region | cr2_region | cr3_region) & core_region ,
|
||||
fill=core)
|
||||
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)
|
||||
|
||||
#Fix control rods initial positions
|
||||
start_pos = 19.2 #cm, geometrical distance between lower bottom and starting point
|
||||
top_pos = 51 * 2.54 # cm, initial position of control rod1 with respect to start post
|
||||
init_pos = 40 * 2.54
|
||||
setattr(cr1_cell, 'translation', [0, 0, start_pos + init_pos])
|
||||
setattr(cr2_cell, 'translation', [-offset, 0, start_pos + top_pos])
|
||||
setattr(cr3_cell, 'translation', [-offset, offset, start_pos + top_pos])
|
||||
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 = 60
|
||||
settings.inactive = 20
|
||||
settings.particles = 500000
|
||||
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)
|
||||
|
||||
#Tallies
|
||||
tally = openmc.Tally(name="heating")
|
||||
if settings.photon_transport:
|
||||
heating_score = 'heating'
|
||||
else:
|
||||
heating_score = 'heating-local'
|
||||
tally.scores.append(heating_score)
|
||||
tallies = openmc.Tallies([tally])
|
||||
|
||||
# Depletion settings
|
||||
model = openmc.model.Model(geometry,mats,settings,tallies)
|
||||
#model.export_to_xml()
|
||||
|
||||
#Plots
|
||||
colors = {salt:'yellow', graphite:'black', inor: 'grey', helium: 'cyan', inconel: 'grey',
|
||||
bush: 'blue', ss316: 'grey', concrete: 'brown', shield: 'red', insulation: 'green',
|
||||
sandwater: 'lightgreen', steel: 'grey'}
|
||||
plots = openmc.Plots()
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xy'
|
||||
plot.width = (150,150)
|
||||
plot.pixels = (1500,1500)
|
||||
plot.origin = (0,0,160)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'xz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (750,1500)
|
||||
plot.origin = (0,-5,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
plot = openmc.Plot()
|
||||
plot.basis = 'yz'
|
||||
plot.width = (150,300)
|
||||
plot.pixels = (500,1000)
|
||||
plot.origin = (5,0,150)
|
||||
plot.color_by = 'material'
|
||||
plot.colors = colors
|
||||
model.plots.append(plot)
|
||||
|
||||
depletion_path = Path(os.getcwd())
|
||||
|
||||
cell = model.geometry.get_cells_by_name('CR1')[0]
|
||||
last_level = get_geom_level_from_res(depletion_path, index=-1)
|
||||
setattr(cell, 'translation', [0, 0, last_level])
|
||||
|
||||
res_path = depletion_path / 'depletion_results.h5'
|
||||
res = openmc.deplete.Results(res_path)
|
||||
|
||||
op = openmc.deplete.CoupledOperator(model,
|
||||
prev_results = res,
|
||||
normalization_mode = "energy-deposition")
|
||||
|
||||
df=pd.read_excel('MSRE_235_233_Power_History_R24E.xlsx')
|
||||
|
||||
#U235
|
||||
timesteps = df['Duration (d)'][:94].values[len(res):]
|
||||
power = df['Power (MWth)'][:94].values[len(res):]*1000000
|
||||
|
||||
#Add one further timestep at the end of every power run
|
||||
# timesteps_ext = []
|
||||
# power_ext = []
|
||||
# for i in range(len(timesteps)):
|
||||
# power_ext.append(power[i])
|
||||
# if power[i] != 0.0:
|
||||
# # duplicate power value
|
||||
# power_ext.append(power[i])
|
||||
# # add timestep - 1 sec
|
||||
# timesteps_ext.append(timesteps[i] - 1/3600/24)
|
||||
# # add 1 sec
|
||||
# timesteps_ext.append(1/3600/24)
|
||||
# else:
|
||||
# timesteps_ext.append(timesteps[i])
|
||||
|
||||
|
||||
|
||||
integrator = openmc.deplete.CECMIntegrator(op, timesteps=timesteps,
|
||||
timestep_units='d', power=power)
|
||||
|
||||
integrator.add_transfer_rate('salt', ['Xe','Kr'], 4.067e-5)
|
||||
integrator.add_transfer_rate('salt', ['Se','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Sb','Te'], 8.777e-3)
|
||||
|
||||
integrator.add_batchwise('trans', axis = 2, cell_id_or_name = 'CR1',
|
||||
bracket = [-2, 5], #cm
|
||||
bracket_limit = [0, top_pos+start_pos],
|
||||
atom_density_limit = 1e-14, #atoms/b-cm
|
||||
tol = 0.1)
|
||||
|
||||
integrator.add_batchwise('refuel', mats_id_or_name = ['salt'],
|
||||
mat_vector = {'U235': 1},
|
||||
bracket = [1e2,1e3], #grams
|
||||
bracket_limit = [0,1e5], #grams
|
||||
tol = 0.01,
|
||||
restart_level = start_pos + init_pos)
|
||||
|
||||
integrator.add_batchwise_wrap('1')
|
||||
|
||||
integrator.integrate()
|
||||
File diff suppressed because one or more lines are too long
|
|
@ -1,17 +1,7 @@
|
|||
# openmc_msre_notebook
|
||||
# OpenMC MSRE_notebooks
|
||||
[](https://www.gnu.org/licenses/gpl-3.0)
|
||||
|
||||
<img src="images/plot_3.png" width="250" height="600"/>
|
||||
<img src="images/lat.png" width="250" height="600"/>
|
||||
|
||||
List of Jupyter Notebooks `Openmc` simulations of the Molten Salt Reactor Experiment (MSRE), operated at ORNL in the 1960s.
|
||||
All scripts are set up using `.h5m` files (surface mesh) of detailed CAD models of the MSRE (designed with `OnShape` CAE tool and available to downlaod [here](https://github.com/openmsr/msre/tree/deplete/step_files)).
|
||||
The surface mesh was created with `Cubit` and is available in the [/h5m](https://github.com/openmsr/openmc_msre_notebooks/tree/main/h5m) folder. Soon a `gmsh` version will also be available.
|
||||
|
||||
## Extra libraries required:
|
||||
- **Pandas**
|
||||
- **Numpy**
|
||||
- **Scypy**
|
||||
- **Seaborn**
|
||||
|
||||
## Depletion msr
|
||||
Some notebooks use a different branch of Openmc: [openmsr/msr](https://github.com/openmsr/openmc/tree/msr_13.2), where msr functionalities have been added.
|
||||
`OpenMC` simulation examples of the Molten Salt Reactor Experiment (MSRE), operated at ORNL in the 1960s.
|
||||
All scripts are set up using the [h5m meshed files](https://github.com/openmsr/msre/tree/master/h5m) obtained with the open source meshing tool [CAD-to-OpenMC](https://github.com/openmsr/CAD_to_OpenMC) from a CAD version of the the MSRE, designed with the CAE tool `OnShape` and available for export here: [onshape msre model]((https://cad.onshape.com/documents/4f04f63bfd4138a61a54b3f8/v/b8c29a0cedda86dfc6948111/)).
|
||||
|
|
|
|||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
Loading…
Add table
Add a link
Reference in a new issue