import openmc
import matplotlib.pyplot as plt
import numpy as np
#Fuel rod regions
central_hole = openmc.ZCylinder(r=0.06)
pallet_out_rad = openmc.ZCylinder(r=0.38)
pallet_in_rad = openmc.ZCylinder(r=0.070)
cladout_rad = openmc.ZCylinder(r=0.455, boundary_type=‘vacuum’)
cladin_rad = openmc.ZCylinder(r=0.3865)
gap = +cladin_rad & -cladout_rad
fuel_max_height = openmc.ZPlane(z0=186.5, boundary_type=‘vacuum’)
fuel_min_height = openmc.ZPlane(z0=-186.50, boundary_type=‘vacuum’)
fuel_height = -fuel_max_height & +fuel_min_height
f_clad_max_height = openmc.ZPlane(z0=206.5, boundary_type=‘vacuum’)
f_clad_height = +fuel_min_height & -f_clad_max_height
mix_fuel_out = openmc.model.HexagonalPrism(edge_length=0.7361215932,orientation=‘x’)
mix_fuel_in = openmc.model.HexagonalPrism(edge_length= 0.6676215932,orientation=‘x’)
mix_fuel_clad_gap = -mix_fuel_out & +mix_fuel_in
rpv_max_height = openmc.ZPlane(z0=559.25)
rpv_min_height = openmc.ZPlane(z0=-559.25)
rpv_height = +rpv_min_height & -rpv_max_height
g_m_channel_out_rad = openmc.ZCylinder(r=0.645)##guiding & measuring channel
g_m_channel_in_rad = openmc.ZCylinder(r=0.545)
ex_r = +central_hole & +fuel_max_height & -f_clad_max_height & -cladin_rad
cell1 = openmc.Cell(region=-pallet_in_rad & +central_hole & fuel_height)
cell2 = openmc.Cell(region=+pallet_out_rad & -cladin_rad & fuel_height)
‘’‘baffle_out_rad = openmc.ZCylinder(r=403.63740285)
baffle_in_rad = openmc.ZCylinder(r=354.63740285)
barrel_out_rad = openmc.ZCylinder(r=418.63740285)
barrel_in_rad = openmc.ZCylinder(r=405.63740285)
rpv_out_rad = openmc.ZCylinder(r=522.13740285)
rpv_in_rad = openmc.ZCylinder(r=482.63740285)’‘’
MATERIALS
Z13_UO2 = openmc.Material(name=“Z13_UO2”)##
Z13_UO2.add_nuclide(“U235”, 0.013)
Z13_UO2.add_nuclide(“U238”, 0.987)
Z13_UO2.add_element(“O”, 2.0)
Z13_UO2.set_density(‘g/cm3’, 10.45)
Z40_UO2 = openmc.Material(name=“Z40_UO2”)##
Z40_UO2.add_nuclide(“U235”, 0.04)
Z40_UO2.add_nuclide(“U238”, 0.96)
Z40_UO2.add_element(“O”, 2.0)
Z40_UO2.set_density(‘g/cm3’, 10.45)
#mix1
Z33Z9_UO2 = openmc.Material(name=“Z33Z9_UO2”)##
Z33Z9_UO2.add_nuclide(“U235”, 0.033)
Z33Z9_UO2.add_nuclide(“U238”, 0.967)
Z33Z9_UO2.add_element(“O”, 2.0)
Z33Z9_UO2.set_density(‘g/cm3’, 10.45)
Z33Z9_UO2_1 = openmc.Material(name=“Z33Z9_UO2_1”)
Z33Z9_UO2_1.add_nuclide(“U235”, 0.024)
Z33Z9_UO2_1.add_nuclide(“U238”, 0.976)
Z33Z9_UO2_1.add_element(“O”, 2.0)
Z33Z9_Gd2O3 = openmc.Material(name=“Z33Z9_Gd2O3”)
Z33Z9_Gd2O3.add_element(“Gd”, 2.0)
Z33Z9_Gd2O3.add_element(“O”, 3.0)
Z33Z9_UO2_Gd2O3 = openmc.Material.mix_materials([Z33Z9_UO2_1, Z33Z9_Gd2O3], [0.92, 0.08], ‘wo’)##
Z33Z9_UO2_Gd2O3.set_density(‘g/cm3’, 10.45)
#mix2
Z44B2_UO2 = openmc.Material(name=“Z44B2_UO2”)##
Z44B2_UO2.add_nuclide(“U235”, 0.044)
Z44B2_UO2.add_nuclide(“U238”, 0.956)
Z44B2_UO2.add_element(“O”, 2.0)
Z44B2_UO2.set_density(‘g/cm3’, 10.45)
Z44B2_UO2_1 = openmc.Material(name=“Z44B2_UO2_1”)
Z44B2_UO2_1.add_nuclide(“U235”, 0.036)
Z44B2_UO2_1.add_nuclide(“U238”, 0.964)
Z44B2_UO2_1.add_element(“O”, 2.0)
Z44B2_Gd2O3 = openmc.Material(name=“Z44B2_Gd2O3”)
Z44B2_Gd2O3.add_element(“Gd”, 2.0)
Z44B2_Gd2O3.add_element(“O”, 3.0)
Z44B2_UO2_Gd2O3 = openmc.Material.mix_materials([Z44B2_UO2_1, Z44B2_Gd2O3], [0.95, 0.05], ‘wo’)##
Z44B2_UO2_Gd2O3.set_density(‘g/cm3’, 10.45)
#mix3
Z33Z2_UO2 = openmc.Material(name=“Z33Z2_UO2”)##
Z33Z2_UO2.add_nuclide(“U235”, 0.033)
Z33Z2_UO2.add_nuclide(“U238”, 0.967)
Z33Z2_UO2.add_element(“O”, 2.0)
Z33Z2_UO2.set_density(‘g/cm3’, 10.45)
Z33Z2_UO2_1 = openmc.Material(name=“Z33Z2_UO2_1”)
Z33Z2_UO2_1.add_nuclide(“U235”, 0.024)
Z33Z2_UO2_1.add_nuclide(“U238”, 0.976)
Z33Z2_UO2_1.add_element(“O”, 2.0)
Z33Z2_Gd2O3 = openmc.Material(name=“Z33Z2_Gd2O3”)
Z33Z2_UO2_1.add_element(“Gd”, 2.0)
Z33Z2_Gd2O3.add_element(“O”, 3.0)
Z33Z2_UO2_Gd2O3 = openmc.Material.mix_materials([Z33Z2_UO2_1, Z33Z2_Gd2O3], [0.92, 0.08], ‘wo’)##
Z33Z2_UO2_Gd2O3.set_density(‘g/cm3’, 10.45)
Z24_UO2 = openmc.Material(name=“Z24_UO2”)##
Z24_UO2.add_nuclide(“U235”, 0.024)
Z24_UO2.add_nuclide(“U238”, 0.976)
Z24_UO2.add_element(“O”, 2.0)
Z24_UO2.set_density(‘g/cm3’, 10.45)
Z24_UO2_1 = openmc.Material(name=“Z24_UO2_1”)
Z24_UO2_1.add_nuclide(“U235”, 0.036)
Z24_UO2_1.add_nuclide(“U238”, 0.964)
Z24_UO2_1.add_element(“O”, 2.0)
Z24_Gd2O3 = openmc.Material(name=“Z24_Gd2O3”)
Z24_Gd2O3.add_element(“Gd”, 2.0)
Z24_Gd2O3.add_element(“O”, 3.0)
Z24_UO2_Gd2O3 = openmc.Material.mix_materials([Z24_UO2_1, Z24_Gd2O3], [0.95, 0.05], ‘wo’)##
Z24_UO2_Gd2O3.set_density(‘g/cm3’, 10.45)
fuel_clad_E110_Zr_Nb = openmc.Material(name=“E110_Zr_Nb”)##clading
fuel_clad_E110_Zr_Nb.add_element(“Zr”, 0.99)
fuel_clad_E110_Zr_Nb.add_element(“Nb”, 0.01)
fuel_clad_E110_Zr_Nb.set_density(‘g/cm3’, 6.515)
gchannel_SS_12H18N10T = openmc.Material(name=“SS_12H18N10T”)##guide_channel_cladding
gchannel_SS_12H18N10T.add_element(“Fe”, 0.68)
gchannel_SS_12H18N10T.add_element(“Cr”, 0.18)
gchannel_SS_12H18N10T.add_element(“Ni”, 0.10)
gchannel_SS_12H18N10T.add_element(“Ti”, 0.01)
gchannel_SS_12H18N10T.add_element(“Mn”, 0.03)
gchannel_SS_12H18N10T.set_density(‘g/cm3’, 7.85)
mchannelE635_Zr_Alloy = openmc.Material(name=“E635_Zr_Alloy”)##measuring_channel_clading
mchannelE635_Zr_Alloy.add_element(“Zr”, 0.935)
mchannelE635_Zr_Alloy.add_element(“Nb”, 0.026)
mchannelE635_Zr_Alloy.add_element(“Sn”, 0.015)
mchannelE635_Zr_Alloy.add_element(“Fe”, 0.024)
mchannelE635_Zr_Alloy.set_density(‘g/cm3’, 6.5)
cbaffel_SS_08H18N10T = openmc.Material(name=“SS_08H18N10T”)##core_baffel_material == core_barel_material
cbaffel_SS_08H18N10T.add_element(“Fe”, 0.70)
cbaffel_SS_08H18N10T.add_element(“Cr”, 0.18)
cbaffel_SS_08H18N10T.add_element(“Ni”, 0.10)
cbaffel_SS_08H18N10T.add_element(“Ti”, 0.01)
cbaffel_SS_08H18N10T.add_element(“Mn”, 0.01)
cbaffel_SS_08H18N10T.set_density(‘g/cm3’, 7.9)
p_vesselSS_15H2NMFA = openmc.Material(name=“SS_15H2NMFA”)##rpv
p_vesselSS_15H2NMFA.add_element(“Fe”, 0.95)
p_vesselSS_15H2NMFA.add_element(“Cr”, 0.015)
p_vesselSS_15H2NMFA.add_element(“Ni”, 0.02)
p_vesselSS_15H2NMFA.add_element(“Mo”, 0.005)
p_vesselSS_15H2NMFA.add_element(“V”, 0.002)
p_vesselSS_15H2NMFA.add_element(“Mn”, 0.008)
p_vesselSS_15H2NMFA.set_density(‘g/cm3’, 7.85)
water = openmc.Material(name=‘Water’)#water
water.add_element(‘H’, 2.0)
water.add_element(‘O’, 1.0)
water.set_density(‘g/cm3’, 1.0)
water.add_s_alpha_beta(‘c_H_in_H2O’)
materials = openmc.Materials([Z13_UO2, Z40_UO2, Z33Z9_UO2, Z33Z9_UO2_Gd2O3, Z44B2_UO2, Z44B2_UO2_Gd2O3, Z33Z2_UO2, Z33Z2_UO2_Gd2O3, Z24_UO2, fuel_clad_E110_Zr_Nb, gchannel_SS_12H18N10T, mchannelE635_Zr_Alloy, cbaffel_SS_08H18N10T, p_vesselSS_15H2NMFA, water, Z24_UO2_Gd2O3])
materials.export_to_xml()
##cladding and central hole cell
central_hole_cell1 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell2 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell3 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell4 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell5 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell6 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell7 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell8 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell9 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
central_hole_cell10 = openmc.Cell(fill=None, region=-central_hole & f_clad_height)
ex_r_cell1 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell2 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell3 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell4 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell5 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell6 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell7 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell8 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell9 = openmc.Cell(fill=None, region=ex_r)
ex_r_cell10 = openmc.Cell(fill=None, region=ex_r)
fuel_clad_cell1 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
fuel_clad_cell2 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
fuel_clad_cell3 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
fuel_clad_cell4 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
fuel_clad_cell5 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
fuel_clad_cell6 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)#
mix_fuel_clad_cell1 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)
mix_fuel_clad_cell2 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)
mix_fuel_clad_cell3 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)
mix_fuel_clad_cell4 = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=gap & f_clad_height)
#mix_fuel_clad_cell = openmc.Cell(fill=fuel_clad_E110_Zr_Nb, region=mix_fuel_clad_gap)
f_outer_cell1 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell2 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell3 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell4 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell5 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell6 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell7 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell8 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell9 = openmc.Cell(fill=water, region=+cladout_rad)
f_outer_cell10 = openmc.Cell(fill=water, region=+cladout_rad)
FA_outer_cell1 = openmc.Cell(fill=water)
FA_outer_cell2 = openmc.Cell(fill=water)
FA_outer_cell3 = openmc.Cell(fill=water)
FA_outer_cell4 = openmc.Cell(fill=water)
FA_outer_cell5 = openmc.Cell(fill=water)
FA_outer_cell6 = openmc.Cell(fill=water)
FA_outer_u1 = openmc.Universe(cells=[FA_outer_cell1])
FA_outer_u2 = openmc.Universe(cells=[FA_outer_cell2])
FA_outer_u3 = openmc.Universe(cells=[FA_outer_cell3])
FA_outer_u4 = openmc.Universe(cells=[FA_outer_cell4])
FA_outer_u5 = openmc.Universe(cells=[FA_outer_cell5])
FA_outer_u6 = openmc.Universe(cells=[FA_outer_cell6])
g_channel_clad = openmc.Cell(fill=gchannel_SS_12H18N10T, region=-g_m_channel_out_rad & +g_m_channel_in_rad & f_clad_height)
m_channel_clad = openmc.Cell(fill=mchannelE635_Zr_Alloy, region=-g_m_channel_out_rad & +g_m_channel_in_rad & f_clad_height)
g_channel_cell = openmc.Cell(region=-g_m_channel_in_rad & f_clad_height)
m_channel_cell = openmc.Cell(region=-g_m_channel_in_rad & f_clad_height)
g_channel_main_cell_u = openmc.Universe(cells=[g_channel_cell, g_channel_clad])
m_channel_main_cell_u = openmc.Universe(cells=[m_channel_cell, m_channel_clad])
Z13_fuel = openmc.Cell(fill=Z13_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z13_fuel_rod_u = openmc.Universe(cells=[Z13_fuel, central_hole_cell1, fuel_clad_cell1, f_outer_cell1, cell1, cell2, ex_r_cell1])
Z33Z9_fuel = openmc.Cell(fill=Z33Z9_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z33Z9_fuel_rod_u = openmc.Universe(cells=[Z33Z9_fuel, central_hole_cell2, fuel_clad_cell2, f_outer_cell2, cell1, cell2, ex_r_cell2])
Z33Z9_mix_fuel = openmc.Cell(fill=Z33Z9_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z33Z9_mix_fuel_rod_u = openmc.Universe(cells=[Z33Z9_mix_fuel, central_hole_cell3, mix_fuel_clad_cell1, f_outer_cell3, cell1, cell2, ex_r_cell3])
Z40_fuel = openmc.Cell(fill=Z40_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z40_fuel_rod_u = openmc.Universe(cells=[Z40_fuel, central_hole_cell4, fuel_clad_cell3, f_outer_cell4, cell1, cell2, ex_r_cell4])
Z44B2_fuel = openmc.Cell(fill=Z44B2_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z44B2_fuel_rod_u = openmc.Universe(cells=[Z44B2_fuel, central_hole_cell5, fuel_clad_cell4, f_outer_cell5, cell1, cell2, ex_r_cell5])
Z44B2_mix_fuel = openmc.Cell(fill=Z24_UO2_Gd2O3, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z44B2_mix_fuel_rod_u = openmc.Universe(cells=[Z44B2_mix_fuel, central_hole_cell6, mix_fuel_clad_cell2, f_outer_cell6, cell1, cell2, ex_r_cell6])
Z24_fuel = openmc.Cell(fill=Z24_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z24_fuel_rod_u = openmc.Universe(cells=[Z24_fuel, central_hole_cell7, fuel_clad_cell5, f_outer_cell7, cell1, cell2, ex_r_cell7])
Z24_mix_fuel = openmc.Cell(fill=Z33Z2_UO2_Gd2O3, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z24_mix_fuel_rod_u = openmc.Universe(cells=[Z24_mix_fuel, central_hole_cell8, mix_fuel_clad_cell3, f_outer_cell8, cell1, cell2, ex_r_cell8])
Z33Z2_fuel = openmc.Cell(fill=Z33Z2_UO2, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z33Z2_fuel_rod_u = openmc.Universe(cells=[Z33Z2_fuel, central_hole_cell9, fuel_clad_cell6, f_outer_cell9, cell1, cell2, ex_r_cell9])
Z33Z2_mix_fuel = openmc.Cell(fill=Z33Z2_UO2_Gd2O3, region=-pallet_out_rad & +pallet_in_rad & fuel_height)
Z33Z2_mix_fuel_rod_u = openmc.Universe(cells=[Z33Z2_mix_fuel, central_hole_cell10, mix_fuel_clad_cell4, f_outer_cell10, cell1, cell2, ex_r_cell10])
outer = openmc.Cell(fill=None)
outer_uni = openmc.Universe(cells=[outer])
#Z13_assembly.pitch = (1.275,)
Z13_assembly = openmc.HexLattice()
Z13_assembly.edge_length = 13.65
Z13_assembly.center = (0, 0)
Z13_assembly.pitch = (1.275,)
Z13_assembly.outer = FA_outer_u1
Z13_assembly.orientation = ‘y’
Z33Z9_assembly = openmc.HexLattice()
Z33Z9_assembly.edge_length = 13.65
Z33Z9_assembly.center = (0, 0)
Z33Z9_assembly.pitch = (1.275,)
Z33Z9_assembly.outer = FA_outer_u2
Z33Z9_assembly.orientation = ‘y’
Z40_assembly = openmc.HexLattice()
Z40_assembly.edge_length = 13.65
Z40_assembly.center = (0, 0)
Z40_assembly.pitch = (1.275,)
Z40_assembly.outer = FA_outer_u3
Z40_assembly.orientation = ‘y’
Z44B2_assembly = openmc.HexLattice()
Z44B2_assembly.edge_length = 13.65
Z44B2_assembly.center = (0, 0)
Z44B2_assembly.pitch = (1.2753,)
Z44B2_assembly.outer = FA_outer_u4
Z44B2_assembly.orientation = ‘y’
Z24_assembly = openmc.HexLattice()
Z24_assembly.edge_length = 13.65
Z24_assembly.center = (0, 0)
Z24_assembly.pitch = (1.3275,)
Z24_assembly.outer = FA_outer_u5
Z24_assembly.orientation = ‘y’
Z33Z2_assembly = openmc.HexLattice()
Z33Z2_assembly.edge_length = 13.65
Z33Z2_assembly.center = (0, 0)
Z33Z2_assembly.pitch = (1.275,)
Z33Z2_assembly.outer = FA_outer_u6
Z33Z2_assembly.orientation = ‘y’
void_boundary = openmc.model.HexagonalPrism(edge_length=13.65, orientation = ‘y’)
void_cell = openmc.Cell(fill=None, region= -void_boundary & f_clad_height)
ring_1 = [Z13_fuel_rod_u]
ring_2 = [Z13_fuel_rod_u]*6
ring_3 = [Z13_fuel_rod_u]*12
ring_4 = [Z13_fuel_rod_u, g_channel_main_cell_u, Z13_fuel_rod_u]*6
ring_5 = [Z13_fuel_rod_u]*11 + [m_channel_main_cell_u] + [Z13_fuel_rod_u]*12
ring_6 = [g_channel_main_cell_u,Z13_fuel_rod_u,Z13_fuel_rod_u,Z13_fuel_rod_u,Z13_fuel_rod_u]*6
ring_7 = [Z13_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z13_fuel_rod_u,Z13_fuel_rod_u,Z13_fuel_rod_u,Z13_fuel_rod_u,Z13_fuel_rod_u,g_channel_main_cell_u]*5 + [Z13_fuel_rod_u,Z13_fuel_rod_u]
ring_8 = [Z13_fuel_rod_u]*42
ring_9 = [Z13_fuel_rod_u]*48
ring_10 = [Z13_fuel_rod_u]*54
ring_11 = [Z13_fuel_rod_u]*60
rring_1 = [Z33Z9_fuel_rod_u]
rring_2 = [Z33Z9_fuel_rod_u]*6
rring_3 = [Z33Z9_fuel_rod_u]*12
rring_4 = [Z33Z9_fuel_rod_u, g_channel_main_cell_u, Z33Z9_fuel_rod_u]*6
rring_5 = [Z33Z9_fuel_rod_u]*5 + [Z33Z9_mix_fuel_rod_u] + [Z33Z9_fuel_rod_u]*5+[m_channel_main_cell_u] + [Z33Z9_fuel_rod_u] + [Z33Z9_mix_fuel_rod_u] + [Z33Z9_fuel_rod_u]*7 + [Z33Z9_mix_fuel_rod_u] + [Z33Z9_fuel_rod_u]*2 #9mix fuel rod position copied from z33z2 assembly
rring_6 = [g_channel_main_cell_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u]*6
rring_7 = [Z33Z9_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z33Z9_fuel_rod_u]*5 +[g_channel_main_cell_u] + [Z33Z9_fuel_rod_u]*5 +[g_channel_main_cell_u] + [Z33Z9_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z33Z9_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z33Z9_fuel_rod_u,]*5 + [g_channel_main_cell_u] + [Z33Z9_fuel_rod_u]*2
rring_8 = [Z33Z9_fuel_rod_u]*42
rring_9 = [Z33Z9_mix_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u,Z33Z9_fuel_rod_u]*6
rring_10 = [Z33Z9_fuel_rod_u]*54
rring_11 = [Z33Z9_fuel_rod_u]*60
rrring_1 = [Z40_fuel_rod_u]
rrring_2 = [Z40_fuel_rod_u]*6
rrring_3 = [Z40_fuel_rod_u]*12
rrring_4 = [Z40_fuel_rod_u,g_channel_main_cell_u,Z40_fuel_rod_u]*6
rrring_5 = [Z40_fuel_rod_u]*11 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*12
rrring_6 = [g_channel_main_cell_u,Z40_fuel_rod_u,Z40_fuel_rod_u,Z40_fuel_rod_u,Z40_fuel_rod_u]*6
rrring_7 = [Z40_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*5 + [m_channel_main_cell_u] + [Z40_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z40_fuel_rod_u]*2
rrring_8 = [Z40_fuel_rod_u]*42
rrring_9 = [Z40_fuel_rod_u]*48
rrring_10 = [Z40_fuel_rod_u]*54
rrring_11 = [Z40_fuel_rod_u]*60
rrrring_1 = [Z44B2_fuel_rod_u]
rrrring_2 = [Z44B2_fuel_rod_u]*6
rrrring_3 = [Z44B2_fuel_rod_u]*12
rrrring_4 = [Z44B2_fuel_rod_u,g_channel_main_cell_u,Z44B2_fuel_rod_u]*6
rrrring_5 = [Z44B2_fuel_rod_u]*5 + [Z44B2_mix_fuel_rod_u] + [Z44B2_fuel_rod_u]*5 + [m_channel_main_cell_u,Z44B2_fuel_rod_u] + [Z44B2_mix_fuel_rod_u] + [Z44B2_fuel_rod_u]*7 + [Z44B2_mix_fuel_rod_u] + [Z44B2_fuel_rod_u]*2
rrrring_6 = [g_channel_main_cell_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u]*6
rrrring_7 = [Z44B2_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z44B2_fuel_rod_u]*2
rrrring_8 = [Z44B2_fuel_rod_u]*42
rrrring_9 = [Z44B2_mix_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u,Z44B2_fuel_rod_u]*6
rrrring_10 = [Z44B2_fuel_rod_u]*54
rrrring_11 = [Z44B2_fuel_rod_u]*60
rrrrring_1 = [Z24_fuel_rod_u]
rrrrring_2 = [Z24_fuel_rod_u]*6
rrrrring_3 = [Z24_fuel_rod_u]*12
rrrrring_4 = [Z24_fuel_rod_u,g_channel_main_cell_u,Z24_fuel_rod_u]*6
rrrrring_5 = [Z24_fuel_rod_u]*11 + [m_channel_main_cell_u] + [Z24_fuel_rod_u]*12
rrrrring_6 = [g_channel_main_cell_u,Z24_fuel_rod_u,Z24_fuel_rod_u,Z24_fuel_rod_u,Z24_fuel_rod_u]*6
rrrrring_7 = [Z24_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*5 + [g_channel_main_cell_u] + [Z24_fuel_rod_u]*2
rrrrring_8 = [Z24_fuel_rod_u]*42
rrrrring_9 = [Z24_mix_fuel_rod_u,Z24_fuel_rod_u,Z24_fuel_rod_u,Z24_fuel_rod_u]*12
rrrrring_10 = [Z24_fuel_rod_u]*54
rrrrring_11 = [Z24_fuel_rod_u]*60
rrrrrring_1 = [Z33Z2_fuel_rod_u]
rrrrrring_2 = [Z33Z2_fuel_rod_u]*6
rrrrrring_3 = [Z33Z2_fuel_rod_u]*12
rrrrrring_4 = [Z33Z2_fuel_rod_u,g_channel_main_cell_u,Z33Z2_fuel_rod_u]*6
rrrrrring_5 = [Z33Z2_fuel_rod_u]*24
rrrrrring_6 = [g_channel_main_cell_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u]*6
rrrrrring_7 = [Z33Z2_fuel_rod_u]*3 + [g_channel_main_cell_u] + [Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,g_channel_main_cell_u]*5 + [Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u]
rrrrrring_8 = [Z33Z2_fuel_rod_u]*42
rrrrrring_9 = [Z33Z2_mix_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u,Z33Z2_fuel_rod_u]*12
rrrrrring_10 = [Z33Z2_fuel_rod_u]*54
rrrrrring_11 = [Z33Z2_fuel_rod_u]*60
Z13_assembly.universes = [
ring_11,
ring_10,
ring_9,
ring_8,
ring_7,
ring_6,
ring_5,
ring_4,
ring_3,
ring_2,
ring_1
]
Z33Z9_assembly.universes = [
rring_11,
rring_10,
rring_9,
rring_8,
rring_7,
rring_6,
rring_5,
rring_4,
rring_3,
rring_2,
rring_1
]
Z40_assembly.universes = [
rrring_11,
rrring_10,
rrring_9,
rrring_8,
rrring_7,
rrring_6,
rrring_5,
rrring_4,
rrring_3,
rrring_2,
rrring_1
]
Z44B2_assembly.universes = [
rrrring_11,
rrrring_10,
rrrring_9,
rrrring_8,
rrrring_7,
rrrring_6,
rrrring_5,
rrrring_4,
rrrring_3,
rrrring_2,
rrrring_1
]
Z24_assembly.universes = [
rrrrring_11,
rrrrring_10,
rrrrring_9,
rrrrring_8,
rrrrring_7,
rrrrring_6,
rrrrring_5,
rrrrring_4,
rrrrring_3,
rrrrring_2,
rrrrring_1
]
Z33Z2_assembly.universes = [
rrrrrring_11,
rrrrrring_10,
rrrrrring_9,
rrrrrring_8,
rrrrrring_7,
rrrrrring_6,
rrrrrring_5,
rrrrrring_4,
rrrrrring_3,
rrrrrring_2,
rrrrrring_1
]
#boundary
h1 = openmc.model.HexagonalPrism(edge_length=13.65, orientation = ‘y’, boundary_type=‘vacuum’)
h2 = openmc.model.HexagonalPrism(edge_length=14., orientation = ‘y’, boundary_type=‘vacuum’)
h = +h1 & -h2
#assembly lattice to cell
Z13_cell = openmc.Cell(fill=Z13_assembly, region=h)
Z33Z9_cell = openmc.Cell(fill=Z33Z9_assembly, region=h)
Z40_cell = openmc.Cell(fill=Z40_assembly, region=h)
Z44B2_cell = openmc.Cell(fill=Z44B2_assembly, region=h)
Z24_cell = openmc.Cell(fill=Z24_assembly, region=h)
Z33Z2_cell = openmc.Cell(fill=Z33Z2_assembly, region=h)
#assembly cell to universe
Z13_assembly_u = openmc.Universe(cells=[Z13_cell])
Z33Z9_assembly_u = openmc.Universe(cells=[Z33Z9_cell])
Z40_assembly_u = openmc.Universe(cells=[Z40_cell])
Z44B2_assembly_u = openmc.Universe(cells=[Z44B2_cell])
Z24_assembly_u = openmc.Universe(cells=[Z24_cell])
Z33Z2_assembly_u = openmc.Universe(cells=[Z33Z2_cell])
void_u = openmc.Universe(cells=[void_cell])
core = openmc.HexLattice()
core.center = (0, 0)
core.pitch = (24.24871131,)#(57.6424935233,)
core.outer = outer_uni
core.orientation = ‘x’
r1 = [Z33Z2_assembly_u]
r2 = [Z13_assembly_u]*6
r3 = [Z24_assembly_u, Z33Z9_assembly_u]*6
r4 = [Z24_assembly_u, Z13_assembly_u, Z13_assembly_u]*6
r5 = [Z13_assembly_u, Z33Z9_assembly_u, Z24_assembly_u, Z33Z9_assembly_u]*6
r6 = [Z24_assembly_u, Z13_assembly_u, Z13_assembly_u, Z13_assembly_u, Z13_assembly_u]*6
r7 = [Z33Z9_assembly_u, Z24_assembly_u, Z33Z2_assembly_u, Z24_assembly_u, Z33Z2_assembly_u, Z24_assembly_u]*6
r8 = [void_u, Z40_assembly_u, Z44B2_assembly_u, Z44B2_assembly_u, Z44B2_assembly_u, Z44B2_assembly_u, Z40_assembly_u]*6
core.universes = [
r8,
r7,
r6,
r5,
r4,
r3,
r2,
r1
]
core_cell = openmc.Cell(fill=core)
core_u = openmc.Universe(cells=[core_cell])#, core_baffle_cell, core_barrel_cell, rpv_cell])
geometry = openmc.Geometry(core_u)
geometry.export_to_xml()
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 5000
Create an initial uniform spatial source distribution over fissionable zones
bounds = [-80, -80, -186.5, 80, 80, 186.5]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.IndependentSource(space=uniform_dist)
Export to “settings.xml”
settings.export_to_xml()
tallies = openmc.Tallies()
mesh = openmc.RegularMesh()
mesh.dimension = [100, 100]
mesh.lower_left = [-100., -100.]
mesh.upper_right = [100., 100.]
Create mesh filter for tally
mesh_filter = openmc.MeshFilter(mesh)
Create mesh tally to score flux and fission rate
tally = openmc.Tally(name=‘flux’)
tally.filters = [mesh_filter]
tally.scores = [‘flux’, ‘fission’]
tallies.append(tally)
tallies.export_to_xml()
openmc.run()
sp = openmc.StatePoint(‘statepoint.100.h5’)
tally = sp.get_tally(scores=[‘flux’])
flux = tally.get_slice(scores=[‘flux’])
fission = tally.get_slice(scores=[‘fission’])
flux.std_dev.shape = (100, 100)
flux.mean.shape = (100, 100)
fission.std_dev.shape = (100, 100)
fission.mean.shape = (100, 100)
fig = plt.subplot(121)
fig.imshow(flux.mean)
fig2 = plt.subplot(122)
fig2.imshow(fission.mean)
relative_error = np.zeros_like(flux.std_dev)
nonzero = flux.mean > 0
relative_error[nonzero] = flux.std_dev[nonzero] / flux.mean[nonzero]
distribution of relative errors
ret = plt.hist(relative_error[nonzero], bins=50)
sp.source[‘E’]
Create log-spaced energy bins from 1 keV to 10 MeV
energy_bins = np.logspace(3,7)
Calculate pdf for source energies
probability, bin_edges = np.histogram(sp.source[‘E’], energy_bins, density=True)
Make sure integrating the PDF gives us unity
print(sum(probability*np.diff(energy_bins)))
Plot source energy PDF
plt.semilogx(energy_bins[:-1], probability*np.diff(energy_bins), drawstyle=‘steps’)
plt.xlabel(‘Energy (eV)’)
plt.ylabel(‘Probability/eV’)
plt.quiver(sp.source[‘r’][‘x’], sp.source[‘r’][‘y’],
sp.source[‘u’][‘x’], sp.source[‘u’][‘y’],
np.log(sp.source[‘E’]), cmap=‘jet’, scale=20.0)
plt.colorbar()
plt.xlim((-0.5,0.5))
plt.ylim((-0.5,0.5))