crp/openmc-designs
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Added some various files

Browse source
This commit is contained in:
Adam Parler 2023-12-16 19:49:46 -08:00
parent 05285cd188
commit 6bc46284d4
6 changed files with 1883 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

274
Depletion/pin_dep.py Normal file
View file

@ -0,0 +1,274 @@
import openmc
import openmc.deplete
import numpy as np
import matplotlib.pyplot as plt
###############################################################################
# Simulation Input File Parameters
###############################################################################
## OpenMC simulation parameters
batches = 300
inactive = 50
particles = 100000
## Depletion simulation parameters
ts = 40*24*60*60 # s
td = 10*24*60*60 # s
time_steps = np.array([ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td])
chain_file = '../chain_endfb71.xml'
power = 17373 # W/cm, for 2D simulations only (use W for 3D)
###############################################################################
# Define materials
###############################################################################
## Instantiate some Materials and register the appropriate Nuclides
zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.49012)
zircaloy.temperature = 600.0
zircaloy.add_nuclide('Sn112', 1.37041E-04, 'wo')
zircaloy.add_nuclide('Sn114', 9.34709E-05, 'wo')
zircaloy.add_nuclide('Sn115', 4.93221E-05, 'wo')
zircaloy.add_nuclide('Sn116', 2.12757E-03, 'wo')
zircaloy.add_nuclide('Sn117', 1.13348E-03, 'wo')
zircaloy.add_nuclide('Sn118', 3.60520E-03, 'wo')
zircaloy.add_nuclide('Sn119', 1.28950E-03, 'wo')
zircaloy.add_nuclide('Sn120', 4.93338E-03, 'wo')
zircaloy.add_nuclide('Sn122', 7.12568E-04, 'wo')
zircaloy.add_nuclide('Sn124', 9.05713E-04, 'wo')
zircaloy.add_nuclide('Fe54', 7.34249E-05, 'wo')
zircaloy.add_nuclide('Fe56', 1.19419E-03, 'wo')
zircaloy.add_nuclide('Fe57', 2.80871E-05, 'wo')
zircaloy.add_nuclide('Fe58', 3.77460E-06, 'wo')
zircaloy.add_nuclide('Cr50', 4.17870E-05, 'wo')
zircaloy.add_nuclide('Cr52', 8.37035E-04, 'wo')
zircaloy.add_nuclide('Cr53', 9.67293E-05, 'wo')
zircaloy.add_nuclide('Cr54', 2.44824E-05, 'wo')
zircaloy.add_nuclide('Ni58', 4.03133E-04, 'wo')
zircaloy.add_nuclide('Ni60', 1.60602E-04, 'wo')
zircaloy.add_nuclide('Ni61', 7.09947E-06, 'wo')
zircaloy.add_nuclide('Ni62', 2.29762E-05, 'wo')
zircaloy.add_nuclide('Ni64', 6.07653E-06, 'wo')
zircaloy.add_nuclide('Zr90', 4.97990E-01, 'wo')
zircaloy.add_nuclide('Zr91', 1.09810E-01, 'wo')
zircaloy.add_nuclide('Zr92', 1.69691E-01, 'wo')
zircaloy.add_nuclide('Zr94', 1.75712E-01, 'wo')
zircaloy.add_nuclide('Zr96', 2.89111E-02, 'wo')
uo2 = openmc.Material(name='Fuel Batch 1')
uo2.set_density('g/cc' ,10.07)
uo2.temperature = 900.0
uo2.add_nuclide('O16' ,1.18530E-01,'wo')
uo2.add_nuclide('U234' ,3.92254E-04,'wo')
uo2.add_nuclide('U235' ,4.36326E-02,'wo')
uo2.add_nuclide('U238' ,8.37445E-01,'wo')
uo2.depletable = True
borated_water = openmc.Material(material_id=4, name='Borated water')
borated_water.set_density('g/cm3', 0.7245)
borated_water.temperature = 600.0
borated_water.add_nuclide('B10',8.40178E-05,'wo')
borated_water.add_nuclide('B11',3.71839E-04,'wo')
borated_water.add_nuclide('H1', 1.11832E-01,'wo')
borated_water.add_nuclide('H2', 3.35290E-05,'wo')
borated_water.add_nuclide('O16',8.85414E-01,'wo' )
borated_water.add_nuclide('O17',2.26386E-03,'wo' )
borated_water.add_s_alpha_beta('c_H_in_H2O')
###############################################################################
# Create geometry
###############################################################################
## Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.41266, name='Fuel OR')
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.474364, name='Clad IR')
left = openmc.XPlane(surface_id=4, x0=-0.632458, name='left')
right = openmc.XPlane(surface_id=5, x0=0.632458, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.632458, name='bottom')
top = openmc.YPlane(surface_id=7, y0=0.632458, name='top')
up = openmc.ZPlane(surface_id=8, z0=50.0, name='up')
down = openmc.ZPlane(surface_id=9, z0=-50.0, name='down')
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
up.boundary_type = 'reflective'
down.boundary_type = 'reflective'
## Instantiate Cells
fuel = openmc.Cell(cell_id=1, name='cell 1')
clad = openmc.Cell(cell_id=3, name='cell 3')
water = openmc.Cell(cell_id=4, name='cell 4')
## Use surface half-spaces to define regions
fuel.region = -fuel_or & -up & +down
clad.region = +fuel_or & -clad_ir & -up & +down
water.region = +clad_ir & +left & -right & +bottom & -top & -up & +down
## Register Materials with Cells
fuel.fill = uo2
clad.fill = zircaloy
water.fill = borated_water
## Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
## Register Cells with Universe
root.add_cells([fuel, clad, water])
## Instantiate a Geometry, register the root Universe
geometry = openmc.Geometry(root)
###############################################################################
# Set volumes of depletable materials
###############################################################################
## Compute cell areas
## area = {}
## area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2
## Set materials volume for depletion. Set to an area for 2D simulations
uo2.volume = 53.49763956
###############################################################################
# Transport calculation settings
###############################################################################
## Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.photon_transport = True
## Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.632458, -0.632458, -50, 0.632458, 0.632458, 50]
# uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.IndependentSource(space=uniform_dist)
## settings_file.electron_treatment = 'ttb'
## settings_file.cutoff = {'energy_photon' : 1000.0}
#entropy_mesh = openmc.RegularMesh()
#entropy_mesh.lower_left = [-0.39218, -0.39218, -50]
#entropy_mesh.upper_right = [0.39218, 0.39218, 50]
#entropy_mesh.dimension = [10, 10, 1]
#settings_file.entropy_mesh = entropy_mesh
###############################################################################
# Initialize and run depletion calculation
###############################################################################
tallies_file = openmc.Tallies()
t1 = openmc.Tally(name='edep')
t1.filters = [openmc.MaterialFilter([uo2,zircaloy,borated_water])]
#t1.filters.append(openmc.ParticleFilter(['neutron']))
#,'photon','electron','positron']))
t1.scores = ['flux','fission', 'nu-fission','heating-local','heating']
tallies_file.append(t1)
tallies_file.export_to_xml()
fission_q={"U232": 193.0442774*1E6,
"U233": 199.7961836*1E6,
"U234": 200.632851*1E6,
"U235": 202.27*1E6,
"U236": 203.4043119*1E6,
"U237": 196.429643*1E6,
"U238": 206.8513817*1E6,
"Np237": 205.3704803*1E6,
"Np238": 208.6993709*1E6,
"Pu238": 209.5400125*1E6,
"Pu239": 208.0185328*1E6,
"Pu240": 208.6125667*1E6,
"Pu241": 211.2377153*1E6,
"Pu242": 212.0721865*1E6,
"Am241": 211.2167986*1E6,
"Am242_m1": 215.1453706*1E6,
"Am243": 212.952779*1E6,
"Cm242": 212.7864913*1E6,
"Cm243": 213.375296*1E6,
"Cm244": 217.9267669*1E6,
"Cm245": 214.6240222*1E6,
"Cm246": 220.179494*1E6,
"Th232": 197.1083894*1E6,
"Pa231": 194.0999425*1E6}
## op = openmc.deplete.Operator(geometry, settings_file, chain_file,fission_q=fission_q)
op = openmc.deplete.Operator(geometry, settings_file, chain_file, energy_mode="energy-deposition")
## Perform simulation using the predictor algorithm
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
integrator.integrate()
## openmc.run(mpi_args=['mpiexec', '-n', '32', '--bind-to', 'core'])
###############################################################################
# Read depletion calculation results
###############################################################################
## Open results file
#results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
## Obtain K_eff as a function of time
#time, keff = results.get_eigenvalue()
## Obtain U235 concentration as a function of time
#time, n_U235 = results.get_atoms('1', 'U235')
## Obtain Xe135 absorption as a function of time
#time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
# t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atom/cm3", time_units="d")
###############################################################################
# Generate plots
###############################################################################
#plt.figure()
#plt.plot(time/(24*60*60), keff, label="K-effective")
#plt.xlabel("Time (days)")
#plt.ylabel("Keff")
#plt.show()
#plt.figure()
#plt.plot(time/(24*60*60), n_U235, label="U 235")
#plt.xlabel("Time (days)")
#plt.ylabel("n U5 (-)")
#plt.show()
#plt.figure()
#plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
#plt.xlabel("Time (days)")
#plt.ylabel("RR (-)")
#plt.show()
#plt.close('all')
#tallies_file = openmc.Tallies()
#energy_filter = openmc.EnergyFilter([0., 20.0e6])
## Instantiate flux Tally in moderator and fuel
#tally = openmc.Tally(name='flux')
#tally.filters = [openmc.CellFilter(fuel)]
#tally.filters.append(energy_filter)
#tally.scores = ['flux']
#tallies_file.append(tally)
## Instantiate reaction rate Tally in moderator
#tally = openmc.Tally(name='fuel rxn rates')
#tally.filters = [openmc.CellFilter(fuel)]
#tally.filters.append(energy_filter)
#tally.scores = ['fission']
#tally.nuclides = ['O16', 'H1']
#tallies_file.append(tally)
#tallies_file.export_to_xml()

119
MSRE/build-it.py Normal file
View file

@ -0,0 +1,119 @@
import openmc
## Materials
## fuel (partially enriched uranium):
## 1200 F pg. 17 MSRE Design and Operations, part III, Nuclear Analysis
temperature = 1200
fuel = openmc.Material(name="MSRE Fuel")
fuel.set_density("g/cm3", 2.146)
fuel.temperature = (1200-32)*5/9
fuel.add_nuclide('Li6', 0.0005, 'wo')
fuel.add_nuclide('Li7', 10.9, 'wo')
fuel.add_nuclide('F19', 66.8, 'wo')
fuel.add_nuclide('Be9', 6.27, 'wo')
fuel.add_element('Zr', 10.92, 'wo')
fuel.add_nuclide('U235', 1.67, 'wo')
fuel.add_nuclide('U238', 3.44, 'wo')
## Moderator graphite
## p. 87 msr operations (robertson) part i
## 1200 F, pg. 17 MSRE Design and Operations, part iii, nuclear analysis
moder = openmc.Material(name="Graphite Moderator")
moder.set_density("g/cc", 1.86)
moder.temperature = (1200-32)*5/9
moder.add_nuclide('B10', 1.592e-5, 'wo')
moder.add_nuclide('B11', 6.408e-5, 'wo')
moder.add_element('V', 0.0009, 'wo')
moder.add_element('S', 0.0005, 'wo')
moder.add_element('C', 99.9852, 'wo')
moder.add_s_alpha_beta(name='C_graphite')
## ignoring the oxygen b/c low content and XS
## Hastelloy Tank
tank = openmc.Material(name="MSRE Tank")
tank.set_density("g/cc", 8.86)
tank.temperature = (1200-32)*5/9
tank.add_nuclide('Si30', 0.00030872, 'wo')
tank.add_nuclide('W186', 0.0014215, 'wo')
tank.add_nuclide('Mn55', 0.008, 'wo')
tank.add_nuclide('W184', 0.001532, 'wo')
tank.add_nuclide('W183', 0.0007155, 'wo')
tank.add_nuclide('W182', 0.001325, 'wo')
tank.add_nuclide('Ni58', 0.4721201092, 'wo')
tank.add_nuclide('Mo92', 0.023744, 'wo')
tank.add_nuclide('Fe58', 0.000141, 'wo')
tank.add_nuclide('Mo94', 0.0148, 'wo')
tank.add_nuclide('Mo95', 0.025472, 'wo')
tank.add_nuclide('Mo96', 0.026688, 'wo')
tank.add_nuclide('Mo97', 0.01528, 'wo')
tank.add_nuclide('Mo98', 0.038608, 'wo')
tank.add_nuclide('Fe56', 0.045877, 'wo')
tank.add_nuclide('Fe57', 0.0010595, 'wo')
tank.add_nuclide('Fe54', 0.0029225, 'wo')
tank.add_nuclide('Mo100', 0.015408, 'wo')
tank.add_nuclide('Si28', 0.00922297, 'wo')
tank.add_nuclide('Si29', 0.00046832, 'wo')
tank.add_nuclide('Cr50', 0.0030415, 'wo')
tank.add_nuclide('Cr52', 0.0586523, 'wo')
tank.add_nuclide('Cr53', 0.0066507, 'wo')
tank.add_nuclide('Cr54', 0.0016555, 'wo')
tank.add_nuclide('Cu65', 0.00107905, 'wo')
tank.add_nuclide('Cu63', 0.00242095, 'wo')
tank.add_nuclide('Ni64', 0.0064191286, 'wo')
tank.add_nuclide('Ni61', 0.0079053205, 'wo')
tank.add_nuclide('Ni60', 0.1818598208, 'wo')
tank.add_nuclide('Ni62', 0.025205621, 'wo')
## Control rod absorber material
## operations report i; p. 102
## 70wtpct gad III oxide, 30wtpct al III oxide
ctrlPois = openmc.Material(name='Control Rod Absorbers')
ctrlPois.set_density("g/cc", 5.873)
ctrlPois.temperature = (1200-32)*5/9
ctrlPois.add_element('Gd', 0.003862068965517241, 'ao')
ctrlPois.add_nuclide('Al27', 0.005884660651235779, 'ao')
ctrlPois.add_nuclide('O16', 0.014620094425129529, 'ao')
air = openmc.Material(name='Air')
air.set_density("g/cc", 0.001225)
air.temperature = (1200-32)*5/9
air.add_nuclide('O16', 0.21, 'ao')
air.add_nuclide('N14', 0.79, 'ao')
surf11 = openmc.Surface(id='11',bound=inf)
surf12 = openmc.Surface(id='12',bound=inf)
surf13 = openmc.Surface(id='13',bound=inf)
surf14 = openmc.Surface(id='14',bound=inf)
cell1 = openmc.Cell(id='111', material=fuel, bound=-11)
cell1 = openmc.Cell(id='112', material=moder, bound=-12)
cell1 = openmc.Cell(id='113', material=tank, bound=-13)
cell1 = openmc.Cell(id='114', material=ctrlPois, bound=-14)
cell1 = openmc.Cell(id='115', material=air, bound=-15)
##### Single Unit Cell
## 4 Bounding planes of stringer
surf1 = openmc.XPlane(x0=-2.54)
surf2 = openmc.XPlane(x0=2.54)
surf3 = openmc.YPlane(y0=-2.54)
surf4 = openmc.YPlane(y0=2.54)
## 4 square cylinders form slot edges
surf5 = openmc.ZCylinder()
## MSR Cluster
msr = openmc.Universe(name='MSR2G-partially-enriched-U-full-core')

306
Rep_Na_3_55.py Normal file
View file

@ -0,0 +1,306 @@
## BISON depletion capability comparison
## Rep_Na_3 assessment case, fuel pellet mesh with 55 radial divisions.
## BISON solves the depletion for regions with the same radius but different volumes.
## The number of radial divisions in the mesh controls the number of the depletion regions
try:
from mpi4py import MPI
rank = MPI.COMM_WORLD.Get_rank()
except ModuleNotFoundError:
rank = 0
import math, sys
import openmc
import openmc.deplete
import numpy as np
from matplotlib import pyplot
######################################################################
# Operational Parameters
######################################################################
fuel_temp = 978.31
gap_temp = 663.69
water_den = 0.7051
ppm_boron = 500
water_temp = 585.25
######################################################################
# Geometric Parameters
######################################################################
Rfo = 0.40956 ## pellet outer radius
Rci = 0.41776 ## clad inner radius
Rco = 0.47736 ## clad outer radius
hpitch = 0.63 ## half-pitch between rods
fheight = 1.3589 ## fuel height
nr = 55 ## number of radial divisions
######################################################################
# Chain Parameters
######################################################################
#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_simple.xml"
#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_casl_pwr.xml"
#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_endfb71_pwr.xml"
chain_path = "/opt/xdata/endfb-vii.1-hdf5/chain_endfb71_pwr.xml"
######################################################################
# Settings
######################################################################
batches = 100
inactive = 20
particles = 50000
seed = np.random.randint(2147483647)*2+1
######################################################################
# Setting the materials
######################################################################
## fuel = 3.0 wt% enriched uranium dioxide
fuel = openmc.Material(material_id=1,name="UO2")
fuel.add_nuclide("U238",0.84223593,'wo')
fuel.add_nuclide("U235",0.03918626,'wo')
fuel.add_nuclide("O16",0.11857781,'wo')
fuel.set_density("g/cc",10.394075)
fuel.temperature = fuel_temp
fuel.volume = (math.pi*Rfo**2)*fheight
mat_list = []
for i in range(nr):
mat = fuel.clone()
mat.volume = (math.pi*((((i+1)/55)*Rfo)**2-(((i)/55)*Rfo)**2))*fheight
mat_list.append(mat)
## gap = helium
gap = openmc.Material(name='He')
gap.set_density('g/cc',0.001598)
gap.add_element('He', 1.0,'wo')
gap.temperature = gap_temp
mat_list.append(gap)
## clad = zircaloy
clad = openmc.Material(name='Zircaloy_4')
clad.set_density('g/cc', 6.56)
clad.add_element('Sn', 0.014 , 'wo')
clad.add_element('Fe', 0.00165, 'wo')
clad.add_element('Cr', 0.001 , 'wo')
clad.add_element('Zr', 0.98335, 'wo')
mat_list.append(clad)
## coolant = borated water
coolant = openmc.Material(name='Borated_water')
coolant.set_density('g/cc', water_den)
coolant.add_nuclide('B10',ppm_boron*0.199E-06,'wo')
coolant.add_nuclide('B11',ppm_boron*0.801E-06,'wo')
coolant.add_element('H',0.1108612,'wo')
coolant.add_element('O',0.8886388,'wo')
coolant.remove_nuclide('O18')
coolant.temperature = water_temp
coolant.add_s_alpha_beta('c_H_in_H2O')
mat_list.append(coolant)
sys.stdout.flush()
try:
MPI.COMM_WORLD.Barrier()
except NameError:
pass
if rank == 0:
print(mat_list)
sys.stdout.flush()
materials = openmc.Materials(mat_list)
materials.cross_sections="/opt/xdata/endfb-vii.1-hdf5/cross_sections.xml"
materials.export_to_xml()
######################################################################
# Building the geometry
######################################################################
## Fuel rings, gap, and cladding
cylinder = []
counter_1 = 0
for i in range(nr):
c = openmc.ZCylinder(surface_id=(i+1), x0=0, y0=0, r=((i+1)/55)*Rfo, name='Fuel '+ str(i+1) + ' OR')
cylinder.append(c)
counter_1 +=1
cylinder_g = openmc.ZCylinder(surface_id=counter_1+1, x0=0, y0=0, r=Rci, name='Gap OR')
cylinder_c = openmc.ZCylinder(surface_id=counter_1+2, x0=0, y0=0, r=Rco, name='Clad OR')
cylinder.append(cylinder_g)
cylinder.append(cylinder_c)
## Box surrounding fuel rod
left = openmc.XPlane(surface_id=counter_1+3, x0=-hpitch, name='left')
right = openmc.XPlane(surface_id=counter_1+4, x0=hpitch, name='right')
bottom = openmc.YPlane(surface_id=counter_1+5, y0=-hpitch, name='bottom')
top = openmc.YPlane(surface_id=counter_1+6, y0=hpitch, name='top')
## Surfaces above/below fuel
fsouth = openmc.ZPlane(surface_id=counter_1+7, z0=-fheight/2.0, name='fsouth')
fnorth = openmc.ZPlane(surface_id=counter_1+8, z0=fheight/2.0, name='fnorth')
## Boundary Conditions
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
fsouth.boundary_type = 'reflective'
fnorth.boundary_type = 'reflective'
cell = []
counter_2 = 0
for i in range(nr):
f = openmc.Cell(cell_id=(i+1), name='fuel '+str(i+1))
cell.append(f)
counter_2 += 1
gap_cell = openmc.Cell(cell_id=counter_2+1, name='Gap cell')
clad_cell = openmc.Cell(cell_id=counter_2+2, name='Clad cell')
water_cell = openmc.Cell(cell_id=counter_2+3, name='Water cell')
cell.append(gap_cell)
cell.append(clad_cell)
cell.append(water_cell)
for i in range(len(cell)-1):
if i==0:
cell[i].region = -cylinder[i] & -fnorth & +fsouth
elif i==(len(cell)-2):
cell[i+1].region = +cylinder[i] & +left & -right & +bottom & -top & -fnorth & +fsouth
cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
else:
cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
count_3 = 0
for i in range(nr):
cell[i].fill = mat_list[i]
count_3 +=1
cell[count_3].fill = gap
cell[count_3+1].fill = clad
cell[count_3+2].fill = coolant
## Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
## Register Cells with Universe
root.add_cells(cell)
## Instantiate a Geometry, register the root Universe, and export to XML
geometry = openmc.Geometry(root)
geometry.export_to_xml()
plot = openmc.Plot.from_geometry(geometry)
plot.pixels = (250, 250)
plot.filename = 'pinplot_55_regions'
plot.width = (hpitch*2, hpitch*2)
plot.pixels = (500, 500)
plot.color_by = 'cell'
plot.colors = {cell[-3]: 'purple',cell[-2]: 'grey',
cell[-1]: 'blue'}
plots = openmc.Plots([plot])
plots.export_to_xml()
openmc.plot_geometry()
###############################################################################
# Exporting to OpenMC settings.xml file
###############################################################################
## Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
settings.seed = seed
settings.verbosity = 7
settings.photon_transport = False
## Cross-section temperature information
settings.temperature = {
'method': 'interpolation',
'tolerance': 100 ,
'multipole': True
}
## Create an initial uniform spatial source distribution over fissionable zones
bounds = [-hpitch, -hpitch, -fheight/2, hpitch, hpitch, fheight/2]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.IndependentSource(space=uniform_dist)
entropy_mesh = openmc.RegularMesh()
entropy_mesh.lower_left = [-hpitch, -hpitch, -1.e50]
entropy_mesh.upper_right = [hpitch, hpitch, 1.e50]
entropy_mesh.dimension = [10, 10, 1]
settings.entropy_mesh = entropy_mesh
settings.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml file
###############################################################################
## Instantiate a tally mesh
mesh = openmc.RegularMesh()
mesh.type = 'regular'
mesh.dimension = [1, 1, 100]
mesh.lower_left = [-hpitch, -hpitch, -fheight/2]
mesh.upper_right = [hpitch, hpitch, fheight/2]
## Instantiate some tally Filters
energy_filter = openmc.EnergyFilter([0., 0.625, 20.e6])
energy_filter2 = openmc.EnergyFilter(np.logspace(-3, 7, num=100))
mesh_filter = openmc.MeshFilter(mesh)
## Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='Spatial Tally 1')
tally.filters = [energy_filter, mesh_filter]
tally.scores = ['flux']#, 'fission', 'nu-fission']
## Instantiate the Tally
tallyspec = openmc.Tally(tally_id=2, name='Energy Spectrum Tally 1')
tallyspec.filters = [energy_filter2]
tallyspec.scores = ['flux']#, 'fission', 'nu-fission']
particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
mat_filter = openmc.MaterialFilter(materials)
tallydep = openmc.Tally(tally_id=3, name='Spatial Tally 2')
tallydep.filters = [mat_filter, particle_filter]
tallydep.scores = ['fission', 'heating', 'heating-local']
## Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([tally, tallyspec,tallydep])
tallies.export_to_xml()
###############################################################################
# Depletion settings
###############################################################################
model = openmc.model.Model(geometry=geometry, materials=materials, settings=settings, tallies=tallies, plots=plots)
## setting the transport operator
operator = openmc.deplete.CoupledOperator(model,chain_path,diff_burnable_mats=False,normalization_mode='fission-q',fission_q={"U235": 202.27e6})
## setting the system linear power [W]
power = [45.50,91.01,136.51,182.02,227.52,273.03,324.45,324.45,273.03,273.03,227.52,182.02,136.51,91.01,45.50]
time_steps = [8640.00,8640.00,8640.00,8640.00,8640.00,8640.00,4268160.00,45792000.00,1728000.00,46915200.00,8640.00,8640.00,8640.00,8640.00,8640.00]
sys.stdout.flush()
try:
MPI.COMM_WORLD.Barrier()
except NameError:
pass
if rank == 0:
print(time_steps)
sys.stdout.flush()
## depleting usin a first-order predictor algorithm
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power, timestep_units = 's')
integrator.integrate()

130
mc-performance/my_attempt.py Normal file
View file

@ -0,0 +1,130 @@
#!/usr/bin/env python3
import openmc
## Monte Carlo Performance Benchmark, version 1.2, July 2011
fuel = openmc.Material(name="Fuel", material_id=1)
fuel.set_density('g/cm3', 10.062)
fuel.add_nuclide("U234", 4.9476e-6)
fuel.add_nuclide("U235", 4.8218e-4)
fuel.add_nuclide("U236", 9.0402e-5)
fuel.add_nuclide("U238", 2.1504e-2)
fuel.add_nuclide("Np237", 7.3733e-6)
fuel.add_nuclide("Pu238", 1.5148e-6)
fuel.add_nuclide("Pu239", 1.3955e-4)
fuel.add_nuclide("Pu240", 3.4405e-5)
fuel.add_nuclide("Pu241", 2.1439e-5)
fuel.add_nuclide("Pu242", 3.7422e-6)
fuel.add_nuclide("Am241", 4.5041e-7)
fuel.add_nuclide("Am242_m1", 9.2301e-9)
fuel.add_nuclide("Am243", 4.7878e-7)
fuel.add_nuclide("Cm242", 1.0485e-7)
fuel.add_nuclide("Cm243", 1.4268e-9)
fuel.add_nuclide("Cm244", 8.8756e-8)
fuel.add_nuclide("Cm245", 3.5285e-9)
fuel.add_nuclide("Mo95", 2.6497e-5)
fuel.add_nuclide("Tc99", 3.2772e-5)
fuel.add_nuclide("Ru101", 3.0742e-5)
fuel.add_nuclide("Ru103", 2.3505e-6)
fuel.add_nuclide("Ag109", 2.0009e-6)
fuel.add_nuclide("Xe135", 1.0801e-8)
fuel.add_nuclide("Cs133", 3.4612e-5)
fuel.add_nuclide("Nd143", 2.6078e-5)
fuel.add_nuclide("Nd145", 1.9898e-5)
fuel.add_nuclide("Sm147", 1.6128e-6)
fuel.add_nuclide("Sm149", 1.1627e-7)
fuel.add_nuclide("Sm150", 7.1727e-6)
fuel.add_nuclide("Sm151", 5.4947e-7)
fuel.add_nuclide("Sm152", 3.0221e-6)
fuel.add_nuclide("Eu153", 2.6209e-6)
fuel.add_nuclide("Gd155", 1.5369e-9)
fuel.add_nuclide("O16", 4.5737e-2)
clad = openmc.Material(name='Cladding', material_id=2)
clad.set_density('g/cm3', 5.77)
clad.add_element("Zr", 1.0)
cold_water = openmc.Material(name='Cold borated water', material_id=3)
cold_water.set_density("atom/b-cm", 0.07416)
cold_water.add_nuclide("H1", 2.0)
cold_water.add_nuclide("O16", 1.0)
cold_water.add_nuclide("B10", 6.490e-4)
cold_water.add_nuclide("B11", 2.689e-3)
cold_water.add_s_alpha_beta('c_H_in_H2O')
hot_water = openmc.Material(name='Hot borated water', material_id=4)
hot_water.set_density("atom/b-cm", 0.06614)
hot_water.add_nuclide("H1", 2.0)
hot_water.add_nuclide("O16", 1.0)
hot_water.add_nuclide("B10", 6.490e-4)
hot_water.add_nuclide("B11", 2.689e-3)
hot_water.add_s_alpha_beta('c_H_in_H2O')
rpv_steel = openmc.Material(name='Reactor pressure vessel steel',
material_id=5)
rpv_steel.set_density("g/cm3", 7.9)
lower_rad_ref = openmc.Material(name='Lower radial reflector',
material_id=6)
lower_rad_ref.set_density("g/cm3", 4.32)
upper_rad_ref = openmc.Material(name='Upper radial reflector /'
'Top plate region', material_id=7)
upper_rad_ref.set_density("g/cm3", 4.28)
bot_plate = openmc.Material(name='Bottom plate region', material_id=8)
bot_plate.set_density("g/cm3", 7.184)
bot_nozzle = openmc.Material(name='Bottom nozzle region', material_id=9)
bot_nozzle.set_density("g/cm3", 2.53)
top_nozzle = openmc.Material(name='Top nozzle region', material_id=10)
top_nozzle.set_density("g/cm3", 1.746)
top_fa = openmc.Material(name='Top of fuel assemblies', material_id=11)
top_fa.set_density("g/cm3", 3.044)
bot_fa = openmc.Material(name='Bottom of fuel assemblies', material_id=12)
bot_fa.set_density("g/cm3", 1.762)
# Define the materials file.
materials = openmc.Materials()
materials.default_xs = '71c'
materials += [fuel, , , , , , , , , , , , bot_fa]
materials.export_to_xml()
# Define surfaces.

401
own_make/build_AP1000.py Normal file
View file

@ -0,0 +1,401 @@
import openmc
from openmc.data import atomic_mass, atomic_weight
MB10 = atomic_mass('B10')
MB11 = atomic_mass('B11')
MU234 = atomic_mass('U234')
MU235 = atomic_mass('U235')
MU238 = atomic_mass('U238')
materials = openmc.Materials()
air = openmc.Material(name='Air')
air.temperature = 300
air.set_density('g/cc', 0.000616)
air.add_element('O', 0.2095, 'ao')
air.add_element('N', 0.7809, 'ao')
air.add_element('Ar', 0.00933, 'ao')
air.add_element('C', 0.00027, 'ao')
## Fuels
'''
fuel_16 = openmc.Material(name='UO2 (1.6%)')
fuel_16.set_density('g/cm3', 10.31341)
fuel_16.temperature = 300
fuel_16.add_nuclide('U234', 3.0131e-06)
fuel_16.add_nuclide('U235', 3.7503e-04)
fuel_16.add_nuclide('U238', 2.2625e-02)
fuel_16.add_nuclide('O16', 4.5897e-02)
fuel_16.add_nuclide('O17', 1.7436e-05)
fuel_16.add_nuclide('O18', 9.2032e-05)
fuel_24 = openmc.Material(name='UO2 (2.4%)')
fuel_24.set_density('g/cm3', 10.29748)
fuel_24.temperature = 300
fuel_24.add_nuclide('U234', 4.4842e-06)
fuel_24.add_nuclide('U235', 5.5814e-04)
fuel_24.add_nuclide('U238', 2.2407e-02)
fuel_24.add_nuclide('O16', 4.5830e-02)
fuel_24.add_nuclide('O17', 1.7411e-05)
fuel_24.add_nuclide('O18', 9.1898e-05)
fuel_31 = openmc.Material(name='UO2 (3.1%)')
fuel_31.set_density('g/cm3', 10.30166)
fuel_31.temperature = 300
fuel_31.add_nuclide('U234', 5.7987e-06)
fuel_31.add_nuclide('U235', 7.2175e-04)
fuel_31.add_nuclide('U238', 2.2253e-02)
fuel_31.add_nuclide('O16', 4.5853e-02)
fuel_31.add_nuclide('O17', 1.7420e-05)
fuel_31.add_nuclide('O18', 9.1942e-05)
fuel_32 = openmc.Material(name='UO2 (3.2%)')
fuel_32.set_density('g/cm3', 10.34115)
fuel_32.temperature = 300
fuel_32.add_nuclide('U234', 5.9959e-06)
fuel_32.add_nuclide('U235', 7.4630e-04)
fuel_32.add_nuclide('U238', 2.2317e-02)
fuel_32.add_nuclide('O16', 4.6029e-02)
fuel_32.add_nuclide('O17', 1.7487e-05)
fuel_32.add_nuclide('O18', 9.2296e-05)
fuel_34 = openmc.Material(name='UO2 (3.4%)')
fuel_34.set_density('g/cm3', 10.35917)
fuel_34.temperature = 300
fuel_34.add_nuclide('U234', 6.4018e-06)
fuel_34.add_nuclide('U235', 7.9681e-04)
fuel_34.add_nuclide('U238', 2.2307e-02)
fuel_34.add_nuclide('O16', 4.6110e-02)
fuel_34.add_nuclide('O17', 1.7517e-05)
fuel_34.add_nuclide('O18', 9.2459e-05)
'''
## Specify enrichments to be calculated
enrs = [0.0161006, 0.0239993, 0.0310221, 0.0319547, 0.0340585]
dens = [10.31341, 10.29748, 10.30166, 10.34115, 10.35917]
# Loop around enrichments
for enr, den in zip(enrs, dens):
# Calculate molar mass of Uranium
enr_25 = enr
enr_24 = 0.008*enr_25
enr_28 = 1.0 - (enr_24 + enr_25)
MU = 1.0/(enr_24/MU234 + enr_25/MU235 + enr_28/MU238)
# Determine molar mass of UO2
MUO2 = MU + 2.0*atomic_weight('O')
# Compute weight percent of U in UO2
wUpUO2 = MU/MUO2
# Calculate Uranium atom fraction
a_U = wUpUO2*MUO2/MU
# Calculate Oxygen atom fraction
a_O = (1.0 - wUpUO2)*MUO2/atomic_weight('O')
# Create material
name = enr*100
mat_name = 'Fuel {0:1.1f}%'.format(name)
mat = openmc.Material(name=mat_name)
mat.temperature = 300
mat.set_density('g/cc', den)
mat.add_element('O', a_O, 'ao')
mat.add_element('U', a_U, 'ao', enrichment=enr*100)
materials += mat
borosilicate_glass = openmc.Material(name='Borosilicate Glass')
borosilicate_glass.set_density('g/cm3', 2.26)
borosilicate_glass.temperature = 300
borosilicate_glass.add_nuclide('Al27', 1.7352e-03)
borosilicate_glass.add_nuclide('B10', 9.6506e-04)
borosilicate_glass.add_nuclide('B11', 3.9189e-03)
borosilicate_glass.add_nuclide('O16', 4.6514e-02)
borosilicate_glass.add_nuclide('O17', 1.7671e-05)
borosilicate_glass.add_nuclide('O18', 9.3268e-05)
borosilicate_glass.add_nuclide('Si28', 1.6926e-02)
borosilicate_glass.add_nuclide('Si29', 8.5944e-04)
borosilicate_glass.add_nuclide('Si30', 5.6654e-04)
ag_in_cd_rods = openmc.Material(name='Ag-In-Cd Control Rod')
ag_in_cd_rods.set_density('g/cm3', 10.16)
ag_in_cd_rods.temperature = 300
ag_in_cd_rods.add_nuclide('Ag107', 2.3523e-02)
ag_in_cd_rods.add_nuclide('Ag109', 2.1854e-02)
ag_in_cd_rods.add_nuclide('Cd106', 3.3882e-05)
ag_in_cd_rods.add_nuclide('Cd108', 2.4166e-05)
ag_in_cd_rods.add_nuclide('Cd110', 3.3936e-04)
ag_in_cd_rods.add_nuclide('Cd111', 3.4821e-04)
ag_in_cd_rods.add_nuclide('Cd112', 6.5611e-04)
ag_in_cd_rods.add_nuclide('Cd113', 3.3275e-04)
ag_in_cd_rods.add_nuclide('Cd114', 7.8252e-04)
ag_in_cd_rods.add_nuclide('Cd116', 2.0443e-04)
ag_in_cd_rods.add_nuclide('In113', 3.4219e-04)
ag_in_cd_rods.add_nuclide('In115', 7.6511e-03)
b4c = openmc.Material(name='B4C Control Rod')
b4c.set_density('g/cm3', 1.76)
b4c.temperature = 300
b4c.add_nuclide('B10', 1.5206e-02)
b4c.add_nuclide('B11', 6.1514e-02)
b4c.add_nuclide('C12', 1.8972e-02)
b4c.add_nuclide('C13', 2.1252e-04)
helium = openmc.Material(name='Helium')
helium.set_density('g/cm3', 0.0015981)
helium.temperature = 300
helium.add_nuclide('He3', 4.8089e-10)
helium.add_nuclide('He4', 2.4044e-04)
inconel_718 = openmc.Material(name='Inconel 718')
inconel_718.set_density('g/cm3', 8.2)
inconel_718.temperature = 300
inconel_718.add_nuclide('Cr50', 7.8239e-04)
inconel_718.add_nuclide('Cr52', 1.5088e-02)
inconel_718.add_nuclide('Cr53', 1.7108e-03)
inconel_718.add_nuclide('Cr54', 4.2586e-04)
inconel_718.add_nuclide('Fe54', 1.4797e-03)
inconel_718.add_nuclide('Fe56', 2.3229e-02)
inconel_718.add_nuclide('Fe57', 5.3645e-04)
inconel_718.add_nuclide('Fe58', 7.1392e-05)
inconel_718.add_nuclide('Mn55', 7.8201e-04)
inconel_718.add_nuclide('Ni58', 2.9320e-02)
inconel_718.add_nuclide('Ni60', 1.1294e-02)
inconel_718.add_nuclide('Ni61', 4.9094e-04)
inconel_718.add_nuclide('Ni62', 1.5653e-03)
inconel_718.add_nuclide('Ni64', 3.9864e-04)
inconel_718.add_nuclide('Si28', 5.6757e-04)
inconel_718.add_nuclide('Si29', 2.8820e-05)
inconel_718.add_nuclide('Si30', 1.8998e-05)
ss304 = openmc.Material(name='Stainless Steel 304')
ss304.set_density('g/cm3', 8.03)
ss304.temperature = 300
ss304.add_nuclide('Cr50', 7.6778e-04)
ss304.add_nuclide('Cr52', 1.4806e-02)
ss304.add_nuclide('Cr53', 1.6789e-03)
ss304.add_nuclide('Cr54', 4.1791e-04)
ss304.add_nuclide('Fe54', 3.4620e-03)
ss304.add_nuclide('Fe56', 5.4345e-02)
ss304.add_nuclide('Fe57', 1.2551e-03)
ss304.add_nuclide('Fe58', 1.6703e-04)
ss304.add_nuclide('Mn55', 1.7604e-03)
ss304.add_nuclide('Ni58', 5.6089e-03)
ss304.add_nuclide('Ni60', 2.1605e-03)
ss304.add_nuclide('Ni61', 9.3917e-05)
ss304.add_nuclide('Ni62', 2.9945e-04)
ss304.add_nuclide('Ni64', 7.6261e-05)
ss304.add_nuclide('Si28', 9.5281e-04)
ss304.add_nuclide('Si29', 4.8381e-05)
ss304.add_nuclide('Si30', 3.1893e-05)
zirc_4 = openmc.Material()
zirc_4.set_density('g/cm3', 6.55)
zirc_4.temperature = 300
zirc_4.add_nuclide('Cr50', 3.2962e-06)
zirc_4.add_nuclide('Cr52', 6.3564e-05)
zirc_4.add_nuclide('Cr53', 7.2076e-06)
zirc_4.add_nuclide('Cr54', 1.7941e-06)
zirc_4.add_nuclide('Fe54', 8.6698e-06)
zirc_4.add_nuclide('Fe56', 1.3610e-04)
zirc_4.add_nuclide('Fe57', 3.1431e-06)
zirc_4.add_nuclide('Fe58', 4.1829e-07)
zirc_4.add_nuclide('O16', 3.0744e-04)
zirc_4.add_nuclide('O17', 1.1680e-07)
zirc_4.add_nuclide('O18', 6.1648e-07)
zirc_4.add_nuclide('Sn112', 4.6735e-06)
zirc_4.add_nuclide('Sn114', 3.1799e-06)
zirc_4.add_nuclide('Sn115', 1.6381e-06)
zirc_4.add_nuclide('Sn116', 7.0055e-05)
zirc_4.add_nuclide('Sn117', 3.7003e-05)
zirc_4.add_nuclide('Sn118', 1.1669e-04)
zirc_4.add_nuclide('Sn119', 4.1387e-05)
zirc_4.add_nuclide('Sn120', 1.5697e-04)
zirc_4.add_nuclide('Sn122', 2.2308e-05)
zirc_4.add_nuclide('Sn124', 2.7897e-05)
zirc_4.add_nuclide('Zr90', 2.1828e-02)
zirc_4.add_nuclide('Zr91', 4.7601e-03)
zirc_4.add_nuclide('Zr92', 7.2759e-03)
zirc_4.add_nuclide('Zr94', 7.3734e-03)
zirc_4.add_nuclide('Zr96', 1.1879e-03)
borated_water = openmc.Material()
borated_water.set_density('g/cm3', 0.740582068)
borated_water.temperature = 300
borated_water.add_nuclide('B10', 7.9714e-06)
borated_water.add_nuclide('B11', 3.2247e-05)
borated_water.add_nuclide('H1', 4.9456e-02)
borated_water.add_nuclide('H2', 7.7035e-06)
borated_water.add_nuclide('O16', 2.4673e-02)
borated_water.add_nuclide('O17', 9.3734e-06)
borated_water.add_nuclide('O18', 4.9474e-05)
plate_borated_water = openmc.Material()
plate_borated_water.set_density('g/cm3', 0.981002532)
plate_borated_water.add_nuclide('B10', 1.0559e-05)
plate_borated_water.add_nuclide('B11', 4.2716e-05)
plate_borated_water.add_nuclide('H1', 6.5512e-02)
plate_borated_water.add_nuclide('H2', 1.0204e-05)
plate_borated_water.add_nuclide('O16', 3.2683e-02)
plate_borated_water.add_nuclide('O17', 1.2416e-05)
plate_borated_water.add_nuclide('O18', 6.5535e-05)
plate_stainless = openmc.Material()
plate_stainless.set_density('g/cm3', 3.68384807)
plate_stainless.temperature = 300
plate_stainless.add_nuclide('Cr50', 3.5223e-04)
plate_stainless.add_nuclide('Cr52', 6.7924e-03)
plate_stainless.add_nuclide('Cr53', 7.7020e-04)
plate_stainless.add_nuclide('Cr54', 1.9172e-04)
plate_stainless.add_nuclide('Fe54', 1.5882e-03)
plate_stainless.add_nuclide('Fe56', 2.4931e-02)
plate_stainless.add_nuclide('Fe57', 5.7578e-04)
plate_stainless.add_nuclide('Fe58', 7.6625e-05)
plate_stainless.add_nuclide('Mn55', 8.0762e-04)
plate_stainless.add_nuclide('Ni58', 2.5731e-03)
plate_stainless.add_nuclide('Ni60', 9.9117e-04)
plate_stainless.add_nuclide('Ni61', 4.3085e-05)
plate_stainless.add_nuclide('Ni62', 1.3738e-04)
plate_stainless.add_nuclide('Ni64', 3.4985e-05)
plate_stainless.add_nuclide('Si28', 4.3711e-04)
plate_stainless.add_nuclide('Si29', 2.2195e-05)
plate_stainless.add_nuclide('Si30', 1.4631e-05)
carbon_steel = openmc.Material()
carbon_steel.set_density('g/cm3', 7.8)
carbon_steel.temperature = 300
carbon_steel.add_nuclide('Al27', 4.3523e-05)
carbon_steel.add_nuclide('B10', 2.5833e-06)
carbon_steel.add_nuclide('B11', 1.0450e-05)
carbon_steel.add_nuclide('C12', 1.0442e-03)
carbon_steel.add_nuclide('C13', 1.1697e-05)
carbon_steel.add_nuclide('Ca40', 1.7043e-05)
carbon_steel.add_nuclide('Ca42', 1.1375e-07)
carbon_steel.add_nuclide('Ca43', 2.3734e-08)
carbon_steel.add_nuclide('Ca44', 3.6673e-07)
carbon_steel.add_nuclide('Ca46', 7.0322e-10)
carbon_steel.add_nuclide('Ca48', 3.2875e-08)
carbon_steel.add_nuclide('Cr50', 1.3738e-05)
carbon_steel.add_nuclide('Cr52', 2.6493e-04)
carbon_steel.add_nuclide('Cr53', 3.0041e-05)
carbon_steel.add_nuclide('Cr54', 7.4778e-06)
carbon_steel.add_nuclide('Cu63', 1.0223e-04)
carbon_steel.add_nuclide('Cu65', 4.5608e-05)
carbon_steel.add_nuclide('Fe54', 4.7437e-03)
carbon_steel.add_nuclide('Fe56', 7.4465e-02)
carbon_steel.add_nuclide('Fe57', 1.7197e-03)
carbon_steel.add_nuclide('Fe58', 2.2886e-04)
carbon_steel.add_nuclide('Mn55', 6.4126e-04)
carbon_steel.add_nuclide('Mo100', 2.9814e-05)
carbon_steel.add_nuclide('Mo92', 4.4822e-05)
carbon_steel.add_nuclide('Mo94', 2.8110e-05)
carbon_steel.add_nuclide('Mo95', 4.8567e-05)
carbon_steel.add_nuclide('Mo96', 5.1015e-05)
carbon_steel.add_nuclide('Mo97', 2.9319e-05)
carbon_steel.add_nuclide('Mo98', 7.4327e-05)
carbon_steel.add_nuclide('Nb93', 5.0559e-06)
carbon_steel.add_nuclide('Ni58', 4.0862e-04)
carbon_steel.add_nuclide('Ni60', 1.5740e-04)
carbon_steel.add_nuclide('Ni61', 6.8420e-06)
carbon_steel.add_nuclide('Ni62', 2.1815e-05)
carbon_steel.add_nuclide('Ni64', 5.5557e-06)
carbon_steel.add_nuclide('P31', 3.7913e-05)
carbon_steel.add_nuclide('S32', 3.4808e-05)
carbon_steel.add_nuclide('S33', 2.7420e-07)
carbon_steel.add_nuclide('S34', 1.5368e-06)
carbon_steel.add_nuclide('S36', 5.3398e-09)
carbon_steel.add_nuclide('Si28', 6.1702e-04)
carbon_steel.add_nuclide('Si29', 3.1330e-05)
carbon_steel.add_nuclide('Si30', 2.0653e-05)
carbon_steel.add_nuclide('Ti46', 1.2144e-06)
carbon_steel.add_nuclide('Ti47', 1.0952e-06)
carbon_steel.add_nuclide('Ti48', 1.0851e-05)
carbon_steel.add_nuclide('Ti49', 7.9634e-07)
carbon_steel.add_nuclide('Ti50', 7.6249e-07)
carbon_steel.add_nuclide('V50', 1.1526e-07)
carbon_steel.add_nuclide('V51', 4.5989e-05)
"""
materials += [fuel_16, fuel_24, fuel_31, fuel_32, fuel_34]
"""
materials += [air, borosilicate_glass, ag_in_cd_rods, b4c, helium,
inconel_718, ss304, zirc_4, borated_water, plate_borated_water,
plate_stainless, carbon_steel]
materials.cross_sections="/opt/xdata/cross_sections.xml"
materials.export_to_xml()
## Fuel Pin
fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(x0=0, y0=0, r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
## Upper Fuel Pin Plenum
spring = openmc.ZCylinder(x0=0, y0=0, r=0.06459, name='Spring')
clad_ir = openmc.ZCylinder(x0=0, y0=0, r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
## Empty Guide Tube above Dashpot
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_or = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide OR')
## Empty Guide Tube at Dashpot
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.50419, name='Guide IR')
guide_or = openmc.ZCylinder(x0=0, y0=0, r=0.54610, name='Guide OR')
## Instrument Tube Pin Geometry
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
## Bare Instrument Thimble Pin
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
## BP above Dashpoint
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.21400, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.23051, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.24130, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.42672, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
## BP Plennum
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.21400, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.23051, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.50419, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.54610, name='Guide IR')
## Control Rod Pin Upper Geometry
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.37338, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
## Control Rod Pin Lower Geometry
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38227, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
## Control Rod Pin Spacer
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.37845, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
## Control Rod Plenum Geometry
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.06459, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')

653
own_make/pin_cell.py Normal file
View file

@ -0,0 +1,653 @@
from numbers import Integral
import numpy as np
import openmc
import openmc.model
def pwr_pin_cell():
"""Create a PWR pin-cell model.
This model is a single fuel pin with 2.4 w/o enriched UO2 corresponding to a
beginning-of-cycle condition and borated water. The specifications are from
the `BEAVRS <http://crpg.mit.edu/research/beavrs>`_ benchmark. Note that the
number of particles/batches is initially set very low for testing purposes.
Returns
-------
model : openmc.model.Model
A PWR pin-cell model
"""
model = openmc.model.Model()
# Define materials.
fuel = openmc.Material(name='UO2 (2.4%)')
fuel.set_density('g/cm3', 10.29769)
fuel.add_nuclide('U234', 4.4843e-6)
fuel.add_nuclide('U235', 5.5815e-4)
fuel.add_nuclide('U238', 2.2408e-2)
fuel.add_nuclide('O16', 4.5829e-2)
clad = openmc.Material(name='Zircaloy')
clad.set_density('g/cm3', 6.55)
clad.add_nuclide('Zr90', 2.1827e-2)
clad.add_nuclide('Zr91', 4.7600e-3)
clad.add_nuclide('Zr92', 7.2758e-3)
clad.add_nuclide('Zr94', 7.3734e-3)
clad.add_nuclide('Zr96', 1.1879e-3)
hot_water = openmc.Material(name='Hot borated water')
hot_water.set_density('g/cm3', 0.740582)
hot_water.add_nuclide('H1', 4.9457e-2)
hot_water.add_nuclide('O16', 2.4672e-2)
hot_water.add_nuclide('B10', 8.0042e-6)
hot_water.add_nuclide('B11', 3.2218e-5)
hot_water.add_s_alpha_beta('c_H_in_H2O')
# Define the materials file.
model.materials = (fuel, clad, hot_water)
# Instantiate ZCylinder surfaces
pitch = 1.26
fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
left = openmc.XPlane(x0=-pitch/2, name='left', boundary_type='reflective')
right = openmc.XPlane(x0=pitch/2, name='right', boundary_type='reflective')
bottom = openmc.YPlane(y0=-pitch/2, name='bottom',
boundary_type='reflective')
top = openmc.YPlane(y0=pitch/2, name='top', boundary_type='reflective')
# Instantiate Cells
fuel_pin = openmc.Cell(name='Fuel', fill=fuel)
cladding = openmc.Cell(name='Cladding', fill=clad)
water = openmc.Cell(name='Water', fill=hot_water)
# Use surface half-spaces to define regions
fuel_pin.region = -fuel_or
cladding.region = +fuel_or & -clad_or
water.region = +clad_or & +left & -right & +bottom & -top
# Create root universe
model.geometry.root_universe = openmc.Universe(0, name='root universe')
model.geometry.root_universe.add_cells([fuel_pin, cladding, water])
model.settings.batches = 10
model.settings.inactive = 5
model.settings.particles = 100
model.settings.source = openmc.Source(space=openmc.stats.Box(
[-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
plot = openmc.Plot.from_geometry(model.geometry)
plot.pixels = (300, 300)
plot.color_by = 'material'
model.plots.append(plot)
return model
def pwr_core():
"""Create a PWR full-core model.
This model is the OECD/NEA Monte Carlo Performance benchmark which is a
grossly simplified pressurized water reactor (PWR) with 241 fuel
assemblies. Note that the number of particles/batches is initially set very
low for testing purposes.
Returns
-------
model : openmc.model.Model
Full-core PWR model
"""
model = openmc.model.Model()
# Define materials.
fuel = openmc.Material(1, name='UOX fuel')
fuel.set_density('g/cm3', 10.062)
fuel.add_nuclide('U234', 4.9476e-6)
fuel.add_nuclide('U235', 4.8218e-4)
fuel.add_nuclide('U238', 2.1504e-2)
fuel.add_nuclide('Xe135', 1.0801e-8)
fuel.add_nuclide('O16', 4.5737e-2)
clad = openmc.Material(2, name='Zircaloy')
clad.set_density('g/cm3', 5.77)
clad.add_nuclide('Zr90', 0.5145)
clad.add_nuclide('Zr91', 0.1122)
clad.add_nuclide('Zr92', 0.1715)
clad.add_nuclide('Zr94', 0.1738)
clad.add_nuclide('Zr96', 0.0280)
cold_water = openmc.Material(3, name='Cold borated water')
cold_water.set_density('atom/b-cm', 0.07416)
cold_water.add_nuclide('H1', 2.0)
cold_water.add_nuclide('O16', 1.0)
cold_water.add_nuclide('B10', 6.490e-4)
cold_water.add_nuclide('B11', 2.689e-3)
cold_water.add_s_alpha_beta('c_H_in_H2O')
hot_water = openmc.Material(4, name='Hot borated water')
hot_water.set_density('atom/b-cm', 0.06614)
hot_water.add_nuclide('H1', 2.0)
hot_water.add_nuclide('O16', 1.0)
hot_water.add_nuclide('B10', 6.490e-4)
hot_water.add_nuclide('B11', 2.689e-3)
hot_water.add_s_alpha_beta('c_H_in_H2O')
rpv_steel = openmc.Material(5, name='Reactor pressure vessel steel')
rpv_steel.set_density('g/cm3', 7.9)
rpv_steel.add_nuclide('Fe54', 0.05437098, 'wo')
rpv_steel.add_nuclide('Fe56', 0.88500663, 'wo')
rpv_steel.add_nuclide('Fe57', 0.0208008, 'wo')
rpv_steel.add_nuclide('Fe58', 0.00282159, 'wo')
rpv_steel.add_nuclide('Ni58', 0.0067198, 'wo')
rpv_steel.add_nuclide('Ni60', 0.0026776, 'wo')
rpv_steel.add_nuclide('Mn55', 0.01, 'wo')
rpv_steel.add_nuclide('Cr52', 0.002092475, 'wo')
rpv_steel.add_nuclide('C0', 0.0025, 'wo')
rpv_steel.add_nuclide('Cu63', 0.0013696, 'wo')
lower_rad_ref = openmc.Material(6, name='Lower radial reflector')
lower_rad_ref.set_density('g/cm3', 4.32)
lower_rad_ref.add_nuclide('H1', 0.0095661, 'wo')
lower_rad_ref.add_nuclide('O16', 0.0759107, 'wo')
lower_rad_ref.add_nuclide('B10', 3.08409e-5, 'wo')
lower_rad_ref.add_nuclide('B11', 1.40499e-4, 'wo')
lower_rad_ref.add_nuclide('Fe54', 0.035620772088, 'wo')
lower_rad_ref.add_nuclide('Fe56', 0.579805982228, 'wo')
lower_rad_ref.add_nuclide('Fe57', 0.01362750048, 'wo')
lower_rad_ref.add_nuclide('Fe58', 0.001848545204, 'wo')
lower_rad_ref.add_nuclide('Ni58', 0.055298376566, 'wo')
lower_rad_ref.add_nuclide('Mn55', 0.0182870, 'wo')
lower_rad_ref.add_nuclide('Cr52', 0.145407678031, 'wo')
lower_rad_ref.add_s_alpha_beta('c_H_in_H2O')
upper_rad_ref = openmc.Material(7, name='Upper radial reflector / Top plate region')
upper_rad_ref.set_density('g/cm3', 4.28)
upper_rad_ref.add_nuclide('H1', 0.0086117, 'wo')
upper_rad_ref.add_nuclide('O16', 0.0683369, 'wo')
upper_rad_ref.add_nuclide('B10', 2.77638e-5, 'wo')
upper_rad_ref.add_nuclide('B11', 1.26481e-4, 'wo')
upper_rad_ref.add_nuclide('Fe54', 0.035953677186, 'wo')
upper_rad_ref.add_nuclide('Fe56', 0.585224740891, 'wo')
upper_rad_ref.add_nuclide('Fe57', 0.01375486056, 'wo')
upper_rad_ref.add_nuclide('Fe58', 0.001865821363, 'wo')
upper_rad_ref.add_nuclide('Ni58', 0.055815129186, 'wo')
upper_rad_ref.add_nuclide('Mn55', 0.0184579, 'wo')
upper_rad_ref.add_nuclide('Cr52', 0.146766614995, 'wo')
upper_rad_ref.add_s_alpha_beta('c_H_in_H2O')
bot_plate = openmc.Material(8, name='Bottom plate region')
bot_plate.set_density('g/cm3', 7.184)
bot_plate.add_nuclide('H1', 0.0011505, 'wo')
bot_plate.add_nuclide('O16', 0.0091296, 'wo')
bot_plate.add_nuclide('B10', 3.70915e-6, 'wo')
bot_plate.add_nuclide('B11', 1.68974e-5, 'wo')
bot_plate.add_nuclide('Fe54', 0.03855611055, 'wo')
bot_plate.add_nuclide('Fe56', 0.627585036425, 'wo')
bot_plate.add_nuclide('Fe57', 0.014750478, 'wo')
bot_plate.add_nuclide('Fe58', 0.002000875025, 'wo')
bot_plate.add_nuclide('Ni58', 0.059855207342, 'wo')
bot_plate.add_nuclide('Mn55', 0.0197940, 'wo')
bot_plate.add_nuclide('Cr52', 0.157390026871, 'wo')
bot_plate.add_s_alpha_beta('c_H_in_H2O')
bot_nozzle = openmc.Material(9, name='Bottom nozzle region')
bot_nozzle.set_density('g/cm3', 2.53)
bot_nozzle.add_nuclide('H1', 0.0245014, 'wo')
bot_nozzle.add_nuclide('O16', 0.1944274, 'wo')
bot_nozzle.add_nuclide('B10', 7.89917e-5, 'wo')
bot_nozzle.add_nuclide('B11', 3.59854e-4, 'wo')
bot_nozzle.add_nuclide('Fe54', 0.030411411144, 'wo')
bot_nozzle.add_nuclide('Fe56', 0.495012237964, 'wo')
bot_nozzle.add_nuclide('Fe57', 0.01163454624, 'wo')
bot_nozzle.add_nuclide('Fe58', 0.001578204652, 'wo')
bot_nozzle.add_nuclide('Ni58', 0.047211231662, 'wo')
bot_nozzle.add_nuclide('Mn55', 0.0156126, 'wo')
bot_nozzle.add_nuclide('Cr52', 0.124142524198, 'wo')
bot_nozzle.add_s_alpha_beta('c_H_in_H2O')
top_nozzle = openmc.Material(10, name='Top nozzle region')
top_nozzle.set_density('g/cm3', 1.746)
top_nozzle.add_nuclide('H1', 0.0358870, 'wo')
top_nozzle.add_nuclide('O16', 0.2847761, 'wo')
top_nozzle.add_nuclide('B10', 1.15699e-4, 'wo')
top_nozzle.add_nuclide('B11', 5.27075e-4, 'wo')
top_nozzle.add_nuclide('Fe54', 0.02644016154, 'wo')
top_nozzle.add_nuclide('Fe56', 0.43037146399, 'wo')
top_nozzle.add_nuclide('Fe57', 0.0101152584, 'wo')
top_nozzle.add_nuclide('Fe58', 0.00137211607, 'wo')
top_nozzle.add_nuclide('Ni58', 0.04104621835, 'wo')
top_nozzle.add_nuclide('Mn55', 0.0135739, 'wo')
top_nozzle.add_nuclide('Cr52', 0.107931450781, 'wo')
top_nozzle.add_s_alpha_beta('c_H_in_H2O')
top_fa = openmc.Material(11, name='Top of fuel assemblies')
top_fa.set_density('g/cm3', 3.044)
top_fa.add_nuclide('H1', 0.0162913, 'wo')
top_fa.add_nuclide('O16', 0.1292776, 'wo')
top_fa.add_nuclide('B10', 5.25228e-5, 'wo')
top_fa.add_nuclide('B11', 2.39272e-4, 'wo')
top_fa.add_nuclide('Zr90', 0.43313403903, 'wo')
top_fa.add_nuclide('Zr91', 0.09549277374, 'wo')
top_fa.add_nuclide('Zr92', 0.14759527104, 'wo')
top_fa.add_nuclide('Zr94', 0.15280552077, 'wo')
top_fa.add_nuclide('Zr96', 0.02511169542, 'wo')
top_fa.add_s_alpha_beta('c_H_in_H2O')
bot_fa = openmc.Material(12, name='Bottom of fuel assemblies')
bot_fa.set_density('g/cm3', 1.762)
bot_fa.add_nuclide('H1', 0.0292856, 'wo')
bot_fa.add_nuclide('O16', 0.2323919, 'wo')
bot_fa.add_nuclide('B10', 9.44159e-5, 'wo')
bot_fa.add_nuclide('B11', 4.30120e-4, 'wo')
bot_fa.add_nuclide('Zr90', 0.3741373658, 'wo')
bot_fa.add_nuclide('Zr91', 0.0824858164, 'wo')
bot_fa.add_nuclide('Zr92', 0.1274914944, 'wo')
bot_fa.add_nuclide('Zr94', 0.1319920622, 'wo')
bot_fa.add_nuclide('Zr96', 0.0216912612, 'wo')
bot_fa.add_s_alpha_beta('c_H_in_H2O')
# Define the materials file.
model.materials = (fuel, clad, cold_water, hot_water, rpv_steel,
lower_rad_ref, upper_rad_ref, bot_plate,
bot_nozzle, top_nozzle, top_fa, bot_fa)
# Define surfaces.
s1 = openmc.ZCylinder(r=0.41, surface_id=1)
s2 = openmc.ZCylinder(r=0.475, surface_id=2)
s3 = openmc.ZCylinder(r=0.56, surface_id=3)
s4 = openmc.ZCylinder(r=0.62, surface_id=4)
s5 = openmc.ZCylinder(r=187.6, surface_id=5)
s6 = openmc.ZCylinder(r=209.0, surface_id=6)
s7 = openmc.ZCylinder(r=229.0, surface_id=7)
s8 = openmc.ZCylinder(r=249.0, surface_id=8, boundary_type='vacuum')
s31 = openmc.ZPlane(z0=-229.0, surface_id=31, boundary_type='vacuum')
s32 = openmc.ZPlane(z0=-199.0, surface_id=32)
s33 = openmc.ZPlane(z0=-193.0, surface_id=33)
s34 = openmc.ZPlane(z0=-183.0, surface_id=34)
s35 = openmc.ZPlane(z0=0.0, surface_id=35)
s36 = openmc.ZPlane(z0=183.0, surface_id=36)
s37 = openmc.ZPlane(z0=203.0, surface_id=37)
s38 = openmc.ZPlane(z0=215.0, surface_id=38)
s39 = openmc.ZPlane(z0=223.0, surface_id=39, boundary_type='vacuum')
# Define pin cells.
fuel_cold = openmc.Universe(name='Fuel pin, cladding, cold water',
universe_id=1)
c21 = openmc.Cell(cell_id=21, fill=fuel, region=-s1)
c22 = openmc.Cell(cell_id=22, fill=clad, region=+s1 & -s2)
c23 = openmc.Cell(cell_id=23, fill=cold_water, region=+s2)
fuel_cold.add_cells((c21, c22, c23))
tube_cold = openmc.Universe(name='Instrumentation guide tube, '
'cold water', universe_id=2)
c24 = openmc.Cell(cell_id=24, fill=cold_water, region=-s3)
c25 = openmc.Cell(cell_id=25, fill=clad, region=+s3 & -s4)
c26 = openmc.Cell(cell_id=26, fill=cold_water, region=+s4)
tube_cold.add_cells((c24, c25, c26))
fuel_hot = openmc.Universe(name='Fuel pin, cladding, hot water',
universe_id=3)
c27 = openmc.Cell(cell_id=27, fill=fuel, region=-s1)
c28 = openmc.Cell(cell_id=28, fill=clad, region=+s1 & -s2)
c29 = openmc.Cell(cell_id=29, fill=hot_water, region=+s2)
fuel_hot.add_cells((c27, c28, c29))
tube_hot = openmc.Universe(name='Instrumentation guide tube, hot water',
universe_id=4)
c30 = openmc.Cell(cell_id=30, fill=hot_water, region=-s3)
c31 = openmc.Cell(cell_id=31, fill=clad, region=+s3 & -s4)
c32 = openmc.Cell(cell_id=32, fill=hot_water, region=+s4)
tube_hot.add_cells((c30, c31, c32))
# Set positions occupied by guide tubes
tube_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8, 11, 14,
2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
tube_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8, 8,
11, 11, 11, 11, 11, 13, 13, 14, 14, 14])
# Define fuel lattices.
l100 = openmc.RectLattice(name='Fuel assembly (lower half)', lattice_id=100)
l100.lower_left = (-10.71, -10.71)
l100.pitch = (1.26, 1.26)
l100.universes = np.tile(fuel_cold, (17, 17))
l100.universes[tube_x, tube_y] = tube_cold
l101 = openmc.RectLattice(name='Fuel assembly (upper half)', lattice_id=101)
l101.lower_left = (-10.71, -10.71)
l101.pitch = (1.26, 1.26)
l101.universes = np.tile(fuel_hot, (17, 17))
l101.universes[tube_x, tube_y] = tube_hot
# Define assemblies.
fa_cw = openmc.Universe(name='Water assembly (cold)', universe_id=5)
c50 = openmc.Cell(cell_id=50, fill=cold_water, region=+s34 & -s35)
fa_cw.add_cell(c50)
fa_hw = openmc.Universe(name='Water assembly (hot)', universe_id=7)
c70 = openmc.Cell(cell_id=70, fill=hot_water, region=+s35 & -s36)
fa_hw.add_cell(c70)
fa_cold = openmc.Universe(name='Fuel assembly (cold)', universe_id=6)
c60 = openmc.Cell(cell_id=60, fill=l100, region=+s34 & -s35)
fa_cold.add_cell(c60)
fa_hot = openmc.Universe(name='Fuel assembly (hot)', universe_id=8)
c80 = openmc.Cell(cell_id=80, fill=l101, region=+s35 & -s36)
fa_hot.add_cell(c80)
# Define core lattices
l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200)
l200.lower_left = (-224.91, -224.91)
l200.pitch = (21.42, 21.42)
l200.universes = [
[fa_cw]*21,
[fa_cw]*21,
[fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
[fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
[fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
[fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
[fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
[fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
[fa_cw]*21,
[fa_cw]*21]
l201 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=201)
l201.lower_left = (-224.91, -224.91)
l201.pitch = (21.42, 21.42)
l201.universes = [
[fa_hw]*21,
[fa_hw]*21,
[fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
[fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
[fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
[fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
[fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
[fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
[fa_hw]*21,
[fa_hw]*21]
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
c1 = openmc.Cell(cell_id=1, fill=l200, region=-s6 & +s34 & -s35)
c2 = openmc.Cell(cell_id=2, fill=l201, region=-s6 & +s35 & -s36)
c3 = openmc.Cell(cell_id=3, fill=bot_plate, region=-s7 & +s31 & -s32)
c4 = openmc.Cell(cell_id=4, fill=bot_nozzle, region=-s5 & +s32 & -s33)
c5 = openmc.Cell(cell_id=5, fill=bot_fa, region=-s5 & +s33 & -s34)
c6 = openmc.Cell(cell_id=6, fill=top_fa, region=-s5 & +s36 & -s37)
c7 = openmc.Cell(cell_id=7, fill=top_nozzle, region=-s5 & +s37 & -s38)
c8 = openmc.Cell(cell_id=8, fill=upper_rad_ref, region=-s7 & +s38 & -s39)
c9 = openmc.Cell(cell_id=9, fill=bot_nozzle, region=+s6 & -s7 & +s32 & -s38)
c10 = openmc.Cell(cell_id=10, fill=rpv_steel, region=+s7 & -s8 & +s31 & -s39)
c11 = openmc.Cell(cell_id=11, fill=lower_rad_ref, region=+s5 & -s6 & +s32 & -s34)
c12 = openmc.Cell(cell_id=12, fill=upper_rad_ref, region=+s5 & -s6 & +s36 & -s38)
root.add_cells((c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12))
# Assign root universe to geometry
model.geometry.root_universe = root
model.settings.batches = 10
model.settings.inactive = 5
model.settings.particles = 100
model.settings.source = openmc.Source(space=openmc.stats.Box(
[-160, -160, -183], [160, 160, 183]))
plot = openmc.Plot()
# plot.basis = 'yz'
plot.type = 'voxel'
plot.origin = (0, 0, 0)
plot.width = (500, 500, 500)
plot.pixels = (250, 250, 250)
# plot.color_by = 'material'
model.plots.append(plot)
return model
def pwr_assembly():
"""Create a PWR assembly model.
This model is a reflected 17x17 fuel assembly from the the `BEAVRS
<http://crpg.mit.edu/research/beavrs>`_ benchmark. The fuel is 2.4 w/o
enriched UO2 corresponding to a beginning-of-cycle condition. Note that the
number of particles/batches is initially set very low for testing purposes.
Returns
-------
model : openmc.model.Model
A PWR assembly model
"""
model = openmc.model.Model()
# Define materials.
fuel = openmc.Material(name='Fuel')
fuel.set_density('g/cm3', 10.29769)
fuel.add_nuclide('U234', 4.4843e-6)
fuel.add_nuclide('U235', 5.5815e-4)
fuel.add_nuclide('U238', 2.2408e-2)
fuel.add_nuclide('O16', 4.5829e-2)
clad = openmc.Material(name='Cladding')
clad.set_density('g/cm3', 6.55)
clad.add_nuclide('Zr90', 2.1827e-2)
clad.add_nuclide('Zr91', 4.7600e-3)
clad.add_nuclide('Zr92', 7.2758e-3)
clad.add_nuclide('Zr94', 7.3734e-3)
clad.add_nuclide('Zr96', 1.1879e-3)
hot_water = openmc.Material(name='Hot borated water')
hot_water.set_density('g/cm3', 0.740582)
hot_water.add_nuclide('H1', 4.9457e-2)
hot_water.add_nuclide('O16', 2.4672e-2)
hot_water.add_nuclide('B10', 8.0042e-6)
hot_water.add_nuclide('B11', 3.2218e-5)
hot_water.add_s_alpha_beta('c_H_in_H2O')
# Define the materials file.
model.materials = (fuel, clad, hot_water)
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
# Create boundary planes to surround the geometry
pitch = 21.42
min_x = openmc.XPlane(x0=-pitch/2, boundary_type='reflective')
max_x = openmc.XPlane(x0=+pitch/2, boundary_type='reflective')
min_y = openmc.YPlane(y0=-pitch/2, boundary_type='reflective')
max_y = openmc.YPlane(y0=+pitch/2, boundary_type='reflective')
# Create a fuel pin universe
fuel_pin_universe = openmc.Universe(name='Fuel Pin')
fuel_cell = openmc.Cell(name='fuel', fill=fuel, region=-fuel_or)
clad_cell = openmc.Cell(name='clad', fill=clad, region=+fuel_or & -clad_or)
hot_water_cell = openmc.Cell(name='hot water', fill=hot_water, region=+clad_or)
fuel_pin_universe.add_cells([fuel_cell, clad_cell, hot_water_cell])
# Create a control rod guide tube universe
guide_tube_universe = openmc.Universe(name='Guide Tube')
gt_inner_cell = openmc.Cell(name='guide tube inner water', fill=hot_water,
region=-fuel_or)
gt_clad_cell = openmc.Cell(name='guide tube clad', fill=clad,
region=+fuel_or & -clad_or)
gt_outer_cell = openmc.Cell(name='guide tube outer water', fill=hot_water,
region=+clad_or)
guide_tube_universe.add_cells([gt_inner_cell, gt_clad_cell, gt_outer_cell])
# Create fuel assembly Lattice
assembly = openmc.RectLattice(name='Fuel Assembly')
assembly.pitch = (pitch/17, pitch/17)
assembly.lower_left = (-pitch/2, -pitch/2)
# Create array indices for guide tube locations in lattice
template_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8,
11, 14, 2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
template_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8,
8, 11, 11, 11, 11, 11, 13, 13, 14, 14, 14])
# Create 17x17 array of universes
assembly.universes = np.tile(fuel_pin_universe, (17, 17))
assembly.universes[template_x, template_y] = guide_tube_universe
# Create root Cell
root_cell = openmc.Cell(name='root cell', fill=assembly)
root_cell.region = +min_x & -max_x & +min_y & -max_y
# Create root Universe
model.geometry.root_universe = openmc.Universe(name='root universe')
model.geometry.root_universe.add_cell(root_cell)
model.settings.batches = 10
model.settings.inactive = 5
model.settings.particles = 100
model.settings.source = openmc.Source(space=openmc.stats.Box(
[-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
plot = openmc.Plot()
plot.basis = 'yz'
plot.origin = (0, 0, 40000)
plot.width = (25, 250)
plot.pixels = (300, 3000)
plot.color_by = 'material'
model.plots.append(plot)
return model
def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5'):
"""Create a 1D slab model.
Parameters
----------
num_regions : int, optional
Number of regions in the problem, each with a unique MGXS dataset.
Defaults to 1.
mat_names : Iterable of str, optional
List of the material names to use; defaults to ['mat_1', 'mat_2',...].
mgxslib_name : str, optional
MGXS Library file to use; defaults to '2g.h5'.
Returns
-------
model : openmc.model.Model
One-group, 1D slab model
"""
openmc.check_type('num_regions', num_regions, Integral)
openmc.check_greater_than('num_regions', num_regions, 0)
if mat_names is not None:
openmc.check_length('mat_names', mat_names, num_regions)
openmc.check_iterable_type('mat_names', mat_names, str)
else:
mat_names = []
for i in range(num_regions):
mat_names.append('mat_' + str(i + 1))
# # Make Materials
materials_file = openmc.Materials()
macros = []
mats = []
for i in range(len(mat_names)):
macros.append(openmc.Macroscopic('mat_' + str(i + 1)))
mats.append(openmc.Material(name=mat_names[i]))
mats[-1].set_density('macro', 1.0)
mats[-1].add_macroscopic(macros[-1])
materials_file += mats
materials_file.cross_sections = mgxslib_name
# # Make Geometry
rad_outer = 929.45
# Set a cell boundary to exist for every material above (exclude the 0)
rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:]
# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
cells = []
surfs = []
surfs.append(openmc.XPlane(x0=0., boundary_type='reflective'))
for r, rad in enumerate(rads):
if r == len(rads) - 1:
surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum'))
else:
surfs.append(openmc.XPlane(x0=rad))
# Instantiate Cells
cells = []
for c in range(len(surfs) - 1):
cells.append(openmc.Cell())
cells[-1].region = (+surfs[c] & -surfs[c + 1])
cells[-1].fill = mats[c]
# Register Cells with Universe
root.add_cells(cells)
# Instantiate a Geometry, register the root Universe, and export to XML
geometry_file = openmc.Geometry(root)
# # Make Settings
# Instantiate a Settings object, set all runtime parameters
settings_file = openmc.Settings()
settings_file.energy_mode = 'multi-group'
settings_file.tabular_legendre = {'enable': False}
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 1000
# Build source distribution
INF = 1000.
bounds = [0., -INF, -INF, rads[0], INF, INF]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
settings_file.source = openmc.Source(space=uniform_dist)
settings_file.output = {'summary': False}
model = openmc.model.Model()
model.geometry = geometry_file
model.materials = materials_file
model.settings = settings_file
model.xs_data = macros
return model
if __name__ == "__main__":
rtn = pwr_core()
rtn.export_to_xml()
openmc.plot_geometry()