Merge pull request #2 from openmsr/open

Open
This commit is contained in:
Aslak Stubsgaard 2022-08-31 14:27:46 +02:00 committed by GitHub
commit ac9365798d
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
17 changed files with 4417323 additions and 3 deletions

View file

@ -1,7 +1,7 @@
# msre
[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
> detailed cad model of the [msre](https://en.wikipedia.org/wiki/Molten-Salt_Reactor_Experiment) (molten salt reactor experiment), operated by oak ridge national laboratory 1965-69.
detailed cad model of the [msre](https://en.wikipedia.org/wiki/Molten-Salt_Reactor_Experiment) (molten salt reactor experiment), operated by oak ridge national laboratory 1965-69.
## msre core
![](core/docs/msre.png)
@ -11,12 +11,35 @@ the work-in-progress cad model can be found [here](https://cad.onshape.com/docum
note that this work and the cad model is under the GNU General Public License v3.0
## prerequisites
### CAD_to_openMC
[CAD_to_openMC](https://github.com/openmsr/CAD_to_openMC) is an open-source package to convert CAD geometry (in the form of '.step' files) into an openmc-readable h5m file
### openmc
these simulations use [openmc](https://docs.openmc.org/en/stable/). automated source installation scripts for linux can be found [here](https://github.com/openmsr/openmc_install_scripts)
## simulation guide
first, clone the repository
```
git clone https://github.com/openmsr/are.git
```
enter the are folder and run the `run.sh` script
```
cd are
bash run.sh
```
## msre heat exchangers
open-access [master's thesis](https://ltu.diva-portal.org/smash/get/diva2:1546993/FULLTEXT01.pdf) produced by Malcolm Akner about simulations of the heat exchangers of the msre, titled:
> Validating results from the Molten Salt Reactor Experiment by use of turbulent CFD simulations
> - A study of a modified U-tube shell-and-tube primary heat exchanger and radiator with molten salts
Validating results from the Molten Salt Reactor Experiment by use of turbulent CFD simulations
- A study of a modified U-tube shell-and-tube primary heat exchanger and radiator with molten salts
### msre primary heat exchanger
![](heatexchanger/docs/phexcadmodel.png)

31
run.sh Executable file
View file

@ -0,0 +1,31 @@
#!/bin/bash
set -ex
PS3='ARE simulations: '
options=("k eigenvalue" "geometry plot" "neutron flux" "photon flux" "quit")
select opt in "${options[@]}"
do
case $opt in
"k eigenvalue")
echo "running k eigenvalue simulation..." &&
bash ./scripts/k.sh
;;
"geometry plot")
echo "plotting geometry..." &&
bash ./scripts/gp.sh
;;
"neutron flux")
echo "generating neutron flux plot..." &&
bash ./scripts/nf.sh
;;
"photon flux")
echo "generating photon flux plot..." &&
bash ./scripts/pf.sh
;;
"quit")
break
;;
*) echo "invalid option $REPLY";;
esac
done

Binary file not shown.

10
scripts/gp.sh Normal file
View file

@ -0,0 +1,10 @@
#!/bin/bash
set -ex
if test -f "./h5m_files/ARE.h5m"; then
python ./scripts/ARE_geometry_plot.py
else
python ./scripts/step_to_h5m.py
python ./scripts/ARE_geometry_plot.py
fi

11
scripts/gs.sh Normal file
View file

@ -0,0 +1,11 @@
#!/bin/bash
set -ex
if test -f "./h5m_files/ARE_gamma.h5m"; then
python ./scripts/ARE_gamma_spectra.py
else
python ./step_to_h5m/step_to_h5m_photon.py
python ./scripts/ARE_gamma_spectra.py
python ./scripts/ARE_gamma_spectra_plotter.py
fi

11
scripts/k.sh Normal file
View file

@ -0,0 +1,11 @@
#!/bin/bash
set -ex
if test -f "./h5m_files/msre.h5m"
then
python ./scripts/msre.py
else
python ./scripts/step_to_h5m.py &&
python ./scripts/msre.py
fi

222
scripts/materials.py Normal file
View file

@ -0,0 +1,222 @@
import sys
sys.path.append('/opt/openmc/')
import openmc
#import neutronsmaterialmaker as nmm
def define_materials():
#salt = openmc.Material(name="salt", temperature = 903)
# salt.add_nuclide('Li7', 10.9566, 'wo')
# salt.add_element('Be', 6.3492, 'wo')
# salt.add_element('Zr', 11.1013, 'wo')
# salt.add_element('Hf', 0.0001, 'wo')
# salt.add_nuclide('U234', 0.0144, 'wo')
# salt.add_nuclide('U235', 1.4093, 'wo')
# salt.add_nuclide('U236', 0.0059, 'wo')
# salt.add_nuclide('U238', 3.0652, 'wo')
# salt.add_element('Fe', 0.0162, 'wo')
# salt.add_element('Cr', 0.0028, 'wo')
# salt.add_element('Ni', 0.0030, 'wo')
# salt.add_element('O', 0.0490, 'wo')
# salt.add_element('F', 67.0270, 'wo')
# salt.set_density('g/cm3', 2.3223)
#
# #fuel salt LiF-BeF2-ZrF4-UF4 (65-29-5-1)
# salt = openmc.Material(name='salt', temperature = 903)
# salt.set_density('g/cm3',2.3275)
# salt.add_nuclide('Li7',65/2)
# salt.add_element('F',65/2+29*2/3+5*4/5+1*4/5)
# salt.add_element('Be',29*1/3)
# salt.add_element('Zr',5*1/5)
# #salt.add_nuclide('U235',1*1/5*0.93)
# #salt.add_nuclide('U238',1*1/5*0.07)
# salt.add_element('U',1*1/5,enrichment=93)
#Salt definition based on report
Li = openmc.Material(name='Li', temperature = 903)
Li.add_nuclide('Li7',1)
Li.set_density('g/cm3',2.3275)
Be = openmc.Material(name='Be', temperature = 903)
Be.add_element('Be',1)
Be.set_density('g/cm3',2.3275)
Zr = openmc.Material(name='Zr', temperature = 903)
Zr.add_element('Zr',1)
Zr.set_density('g/cm3',2.3275)
Hf = openmc.Material(name='Hf', temperature = 903)
Hf.add_element('Hf',1)
Hf.set_density('g/cm3',2.3275)
U234 = openmc.Material(name='U234', temperature = 903)
U234.add_nuclide('U234',1)
U234.set_density('g/cm3',2.3275)
U235 = openmc.Material(name='U235', temperature = 903)
U235.add_nuclide('U235',1)
U235.set_density('g/cm3',2.3275)
U236 = openmc.Material(name='U236', temperature = 903)
U236.add_nuclide('U236',1)
U236.set_density('g/cm3',2.3275)
U238 = openmc.Material(name='U238', temperature = 903)
U238.add_nuclide('U238',1)
U238.set_density('g/cm3',2.3275)
Fe = openmc.Material(name='Fe', temperature = 903)
Fe.add_element('Fe',1)
Fe.set_density('g/cm3',2.3275)
Cr = openmc.Material(name='Cr', temperature = 903)
Cr.add_element('Cr',1)
Cr.set_density('g/cm3',2.3275)
Ni = openmc.Material(name='Ni', temperature = 903)
Ni.add_element('Ni',1)
Ni.set_density('g/cm3',2.3275)
O = openmc.Material(name='O', temperature = 903)
O.add_element('O',1)
O.set_density('g/cm3',2.3275)
F = openmc.Material(name='F', temperature = 903)
F.add_element('F',1)
F.set_density('g/cm3',2.3275)
salt = openmc.Material.mix_materials(
[Li, Be, Zr, Hf, U234, U235, U236, U238, Fe, Cr, Ni, O, F],
[10.9566/100, 6.3492/100, 11.1013/100, 0.0001/100, 0.0144/100, 1.4093/100, 0.0059/100, 3.0652/100, 0.0162/100, 0.0028/100, 0.0030/100, 0.0490/100, 67.0270/100],
'wo')
salt.name="salt"
#moderator blocks170
graphite = openmc.Material(name='graphite',temperature=903)
graphite.set_density('g/cm3',1.86)
graphite.add_element('C',1)
graphite.add_s_alpha_beta('c_Graphite')
#inor
Mo = openmc.Material(name='Mo', temperature = 903)
Mo.add_element('Mo',1)
Mo.set_density('g/cm3',2.3275)
C = openmc.Material(name='C', temperature = 903)
C.add_element('Mo',1)
C.set_density('g/cm3',2.3275)
Al = openmc.Material(name='Al', temperature = 903)
Al.add_element('Al',1)
Al.set_density('g/cm3',2.3275)
Ti = openmc.Material(name='Ti', temperature = 903)
Ti.add_element('Ti',1)
Ti.set_density('g/cm3',2.3275)
S = openmc.Material(name='S', temperature = 903)
S.add_element('S',1)
S.set_density('g/cm3',2.3275)
Mn = openmc.Material(name='Mn', temperature = 903)
Mn.add_element('Mn',1)
Mn.set_density('g/cm3',2.3275)
Si = openmc.Material(name='Si', temperature = 903)
Si.add_element('Si',1)
Si.set_density('g/cm3',2.3275)
Cu = openmc.Material(name='Cu', temperature = 903)
Cu.add_element('Cu',1)
Cu.set_density('g/cm3',2.3275)
B = openmc.Material(name='B', temperature = 903)
B.add_element('B',1)
B.set_density('g/cm3',2.3275)
W = openmc.Material(name='W', temperature = 903)
W.add_element('W',1)
W.set_density('g/cm3',2.3275)
P = openmc.Material(name='P', temperature = 903)
P.add_element('P',1)
P.set_density('g/cm3',2.3275)
Co = openmc.Material(name='Co', temperature = 903)
Co.add_element('Co',1)
Co.set_density('g/cm3',2.3275)
inor = openmc.Material.mix_materials(
[Ni, Mo, Cr, Fe, C, Al, Ti, S, Mn, Si, Cu, B, W, P, Co],
[68/100, 17/100, 7/100, 5/100, 0.06/100, 0.2/100, 0.2/100, 0.02/100, 0.8/100, 0.8/100, 0.3/100, 0.01/100, 0.4/100, 0.015/100, 0.195/100],
'wo')
inor.name="inor"
inor.set_density('g/cm3',8.7745)
# inor = openmc.Material(name='inor',temperature=903)
# 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 = 903)
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)
#Control rods bushing posion material
Gd2O3 = openmc.Material()
Gd2O3.add_element('Gd',2)
Gd2O3.add_element('O',3)
Gd2O3.set_density('g/cm3',5.873)
Al2O3 = openmc.Material()
Al2O3.add_element('Al',2)
Al2O3.add_element('O',3)
Al2O3.set_density('g/cm3',5.873)
bush = openmc.Material.mix_materials([Gd2O3,Al2O3],[0.7,0.3],'wo')
bush.name='bush'
#Concrete block
#concrete=nmm.Material.from_library('Concrete, Regular')
#concrete=concrete.openmc_material
#concrete.name="concrete"
#Baryte block
baryte = openmc.Material(name="baryte", temperature = 903)
baryte.add_element('Ba',1)
baryte.add_element('S',1)
baryte.add_element('O',4)
baryte.set_density('g/cm3',4.5)
#Thermal shielding
water = openmc.Material()
water.add_element('H',2)
water.add_element('O',1)
water.set_density('g/cm3',0.997)
#steel = nmm.Material.from_library('Steel, Carbon')
#steel=steel.openmc_material
#shield = openmc.Material.mix_materials([water,steel],[0.5,0.5],'wo')
#shield.name='waterSteel'
detector = openmc.Material(name = 'water')
detector.add_element('O',76.2,percent_type='wo')
detector.add_element('C',11.1,percent_type='wo')
detector.add_element('H',10.1,percent_type='wo')
detector.add_element('N',2.6,percent_type='wo')
detector.set_density('g/cm3',1.0)
mats = openmc.Materials([salt,graphite,inor,helium,inconel,bush])#,concrete,baryte,shield,detector])
mats.export_to_xml()
return mats

36
scripts/msre.py Normal file
View file

@ -0,0 +1,36 @@
import matplotlib.pyplot as plt
from materials import *
import openmc
###############################################################################
#k eigenvalue simulation on ARE (all safety rods fully inserted)
###############################################################################
h5m_filepath = 'h5m_files/msre.h5m'
#Materials
define_materials()
#geometry
graveyard=openmc.Sphere(r=10000,boundary_type='vacuum')
cad_univ = openmc.DAGMCUniverse(filename=h5m_filepath,auto_geom_ids=True,universe_id=996 )
cad_cell = openmc.Cell(cell_id=997 , region= -graveyard, fill= cad_univ)
root = openmc.Universe(universe_id=998)
root.add_cells([cad_cell])
geometry = openmc.Geometry(root)
geometry.export_to_xml()
settings = openmc.Settings()
settings.temperature = {'method':'interpolation'}
settings.batches = 100
settings.inactive = 10
settings.particles = 10000
settings.max_lost_particles = 10000
source_area = openmc.stats.Box([-200., -200., -200.],[ 200., 200., 200.],only_fissionable = True)
settings.source = openmc.Source(space=source_area)
settings.export_to_xml()
openmc.run()

View file

@ -0,0 +1,57 @@
import matplotlib
import openmc
from materials import *
###############################################################################
#generate geometry plot of are (all safety rods fully inserted)
###############################################################################
#geometry
h5m_filepath = '../h5m_files/msre_simp_cubit/msre_simp_1e-2.h5m'
graveyard = openmc.Sphere(r=10000,boundary_type='vacuum')
cad_univ = openmc.DAGMCUniverse(filename=h5m_filepath,auto_geom_ids=True)
cad_cell = openmc.Cell(region=-graveyard,fill=cad_univ)
root=openmc.Universe()
root.add_cells([cad_cell])
geometry=openmc.Geometry(root)
geometry.export_to_xml()
# materials
mats = define_materials()
#mats = openmc.Materials([salt,BeO,inconel,insulation,coolant,helium,stainless,boron])
mats.export_to_xml()
#plotting geometry
plots = openmc.Plots()
x_width = 300
y_width = 300
#xy plot
p1 = openmc.Plot()
p1.width = (x_width,y_width)
p1.origin = (0,0,100)
p1.pixels = (1000, 1000)
p1.color_by = 'material'
#xz plot (split plane)
p2 = openmc.Plot()
p2.basis = 'xz'
p2.origin = (0,0,130)
p2.width = (x_width,y_width)
p2.pixels = (1000, 1000)
p2.color_by = 'material'
p3 = openmc.Plot()
p3.basis = 'yz'
p3.origin = (0,0,130)
p3.width = (x_width,y_width)
p3.pixels = (1000, 1000)
p3.color_by = 'material'
plots.append(p1)
plots.append(p2)
plots.append(p3)
plots.export_to_xml()
openmc.plot_geometry()

View file

@ -0,0 +1,68 @@
import matplotlib.pyplot as plt
import openmc
from materials import *
###############################################################################
#create .png file of neutron flux (all safety rods inserted)
###############################################################################
#Geometry
h5m_filepath = '../h5m_files/msre_simp_cubit/msre_simp_1e-2.h5m'
graveyard = openmc.Sphere(r=10000,boundary_type='vacuum')
cad_univ = openmc.DAGMCUniverse(filename=h5m_filepath,auto_geom_ids=True)
cad_cell = openmc.Cell(region=-graveyard,fill=cad_univ)
root=openmc.Universe()
root.add_cells([cad_cell])
geometry=openmc.Geometry(root)
geometry.export_to_xml()
#materials
mats=define_materials()
#mats = openmc.Materials([salt,BeO,inconel,insulation,coolant,helium,stainless,boron,blanket,shield])
mats.export_to_xml()
settings = openmc.Settings()
settings.temperature = {'method':'interpolation'}
settings.batches = 100
settings.inactive = 10
settings.particles = 5000
settings.max_lost_particles = 10000
source_area = openmc.stats.Box([-200., -200., -200.],[ 200., 200., 200.],only_fissionable = True)
settings.source = openmc.Source(space=source_area)
settings.export_to_xml()
#tallies
tallies = openmc.Tallies()
mesh = openmc.RegularMesh()
mesh.dimension = [1000,1000]
mesh.lower_left = [-300,-300]
mesh.upper_right = [300,300]
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally(name='flux')
tally.filters = [mesh_filter]
tally.scores = ['flux','fission']
tallies.append(tally)
tallies.export_to_xml()
model = openmc.model.Model(geometry, mats, settings, tallies)
sp_filename = model.run()
sp = openmc.StatePoint(sp_filename)
s_tally = sp.get_tally(scores=['flux','fission'])
flux = s_tally.get_slice(scores=['flux'])
fission = s_tally.get_slice(scores=['fission'])
flux.std_dev.shape = (1000,1000)
flux.mean.shape = (1000,1000)
fission.std_dev.shape = (1000,1000)
fission.mean.shape = (1000,1000)
fig = plt.subplot(121)
fig.axis([350,650,350,650])
fig.pixels = (2000,2000)
fig.imshow(flux.mean)
plt.savefig('neutron_flux', dpi=2000)

View file

@ -0,0 +1,74 @@
import matplotlib.pyplot as plt
import openmc
from materials import *
###############################################################################
#create .png file of photon flux (all safety rods inserted)
###############################################################################
#Geometry
h5m_filepath = 'h5m_files/ARE.h5m'
graveyard = openmc.Sphere(r=10000,boundary_type='vacuum')
cad_univ = openmc.DAGMCUniverse(filename=h5m_filepath,auto_geom_ids=True)
cad_cell = openmc.Cell(region=-graveyard,fill=cad_univ)
root=openmc.Universe()
root.add_cells([cad_cell])
geometry=openmc.Geometry(root)
geometry.export_to_xml()
#materials
mats = openmc.Materials([salt,BeO,inconel,insulation,coolant,helium,stainless,boron,blanket,shield])
mats.export_to_xml()
settings = openmc.Settings()
settings.temperature = {'method':'interpolation'}
settings.batches = 100
settings.inactive = 10
settings.particles = 5000
settings.photon_transport = True
source_area = openmc.stats.Box([-200., -200., -200.],[ 200., 200., 200.],only_fissionable = True)
settings.source = openmc.Source(space=source_area)
settings.export_to_xml()
#tallies
tallies = openmc.Tallies()
# sets up filters for the tallies
photon_particle_filter = openmc.ParticleFilter(['photon']) # note the use of photons here
energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175']
energy_filter = openmc.EnergyFilter(energy_bins)
mesh = openmc.RegularMesh()
mesh.dimension = [1000,1000]
mesh.lower_left = [-300,-300]
mesh.upper_right = [300,300]
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally(name='flux')
tally.filters = [mesh_filter,photon_particle_filter]
tally.scores = ['flux']
tallies.append(tally)
tallies.export_to_xml()
# combine all the required parts to make a model
model = openmc.model.Model(geom, mats, settings, tallies)
# remove old files and runs OpenMC
results_filename = model.run()
sp = openmc.StatePoint('statepoint.100.h5')
s_tally = sp.get_tally(scores=['flux'])
flux = s_tally.get_slice(scores=['flux'])
flux.std_dev.shape = (1000,1000)
flux.mean.shape = (1000,1000)
fig = plt.subplot(121)
fig.axis([350,650,350,650])
fig.pixels = (2000,2000)
fig.imshow(flux.mean)
plt.savefig('photon_flux', dpi=2000)

10
scripts/nf.sh Normal file
View file

@ -0,0 +1,10 @@
#!/bin/bash
set -ex
if test -f "./h5m_files/ARE.h5m"; then
python ./scripts/ARE_neutron_flux.py
else
python ./scripts/step_to_h5m.py
python ./scripts/ARE_neutron_flux.py
fi

10
scripts/pf.sh Normal file
View file

@ -0,0 +1,10 @@
#!/bin/bash
set -ex
if test -f "./h5m_files/ARE.h5m"; then
python ./scripts/ARE_photon_flux.py
else
python ./scripts/step_to_h5m.py
python ./scripts/ARE_photon_flux.py
fi

29
scripts/step_to_h5m.py Normal file
View file

@ -0,0 +1,29 @@
###############################################################################
# Converting step files to h5m file to be read by openmc
###############################################################################
import numpy as np
import os
import sys
sys.path.append('../CAD_to_openMC/src/')
import CAD_to_OpenMC.assembly as ab
###############################################################################
# inputs
step_filepath = "./step_files/msre.step"
h5m_out_filepath = os.getcwd() + '/h5m_files/msre.h5m'
# mesher config
ab.mesher_config['min_mesh_size'] = 0.2
ab.mesher_config['mesh_algorithm'] = 4
ab.mesher_config['threads'] = 6
#ab.mesher_config['max_mesh_size'] = 1000
#ab.mesher_config['curve_samples'] = 100
#ab.mesher_config['vetoed'] = [474,1157,1243,1341,1537]
#ab.mesher_config['angular_tolerance'] = 0.09
#ab.mesher_config['tolerance'] = 0.001
# output
a=ab.Assembly()
a.stp_files=[step_filepath]
a.import_stp_files()
a.solids_to_h5m(backend='gmsh',h5m_filename=h5m_out_filepath)

1345087
step_files/msre.step Normal file

File diff suppressed because it is too large Load diff

1539188
step_files/msre_no_cr.step Normal file

File diff suppressed because it is too large Load diff

1532453
step_files/msre_simp.step Normal file

File diff suppressed because it is too large Load diff