Hello everyone!
i’m running a carbide 3600 program according to benchmark using openmc. when I’m looking for a keff score, the grade generated in the 12th month never gets a score of 1. I’ve also used enrichment up to 20% but it never reaches 1 in the 12th month. Can you help find errors in the program I made? 
program:
import openmc
import math
import matplotlib.pyplot as plt
from matplotlib import pyplot
import openmc.deplete
import pandas as pd
import numpy as np
#Definisi Material#
####Inisialisasi Material/memasukkan material####
#SET_TEMPERATURE
fuel_temperature = (987) #K #987
structure_temperature = (470) #K #470
#0 - 20.11 cm from bottom active core
inner_1 = openmc.Material(name = ‘0 - 20.11 cm from bottom active core’)
inner_1.volume = np.pi*(0.34700.3954)20.1146930
inner_1.temperature = fuel_temperature
inner_1.add_nuclide(‘C12’, 2.6385E-02)
inner_1.add_nuclide(‘U234’, 1.1919E-06)
inner_1.add_nuclide(‘U235’, 4.6799E-05)
inner_1.add_nuclide(‘U236’, 2.8136E-06)
inner_1.add_nuclide(‘U238’, 2.2303E-02)
inner_1.add_nuclide(‘Np237’, 2.5168E-06)
inner_1.add_nuclide(‘Np239’, 3.5946E-06)
inner_1.add_nuclide(‘Pu238’, 1.0646E-04)
inner_1.add_nuclide(‘Pu239’, 1.8140E-03)
inner_1.add_nuclide(‘Pu240’, 1.0229E-03)
inner_1.add_nuclide(‘Pu241’, 2.4210E-04)
inner_1.add_nuclide(‘Pu242’, 3.4136E-04)
inner_1.add_nuclide(‘Am241’, 3.7572E-05)
inner_1.add_nuclide(‘Am242’, 1.0156E-08)
inner_1.add_nuclide(‘Am242’, 6.9199E-07)
inner_1.add_nuclide(‘Am243’, 1.4205E-05)
inner_1.add_nuclide(‘Cm242’, 1.0156E-08)
inner_1.add_nuclide(‘Cm243’, 1.6209E-06)
inner_1.add_nuclide(‘Cm244’, 4.3067E-08)
inner_1.add_nuclide(‘Cm245’, 1.0655E-06)
inner_1.add_nuclide(‘Cm246’, 2.7855E-08)
inner_1.add_element(‘Mo’, 1.3400E-03)
inner_1.set_density(‘sum’)
#20.11 cm - 40.22 cm from bottom active core
inner_2 = openmc.Material(name = ‘20.11 - 40.22 cm from bottom active core’)
inner_2.volume = np.pi*(0.34700.3954)20.1146930
inner_2.temperature = fuel_temperature
inner_2.add_nuclide(‘C12’, 2.6385E-02)
inner_2.add_nuclide(‘U234’, 1.1668E-06)
inner_2.add_nuclide(‘U235’, 4.4631E-05)
inner_2.add_nuclide(‘U236’, 3.1423E-06)
inner_2.add_nuclide(‘U238’, 2.2143E-02)
inner_2.add_nuclide(‘Np237’, 3.4681E-06)
inner_2.add_nuclide(‘Np239’, 4.4531E-06)
inner_2.add_nuclide(‘Pu238’, 1.0306E-04)
inner_2.add_nuclide(‘Pu239’, 1.8594E-03)
inner_2.add_nuclide(‘Pu240’, 1.0161E-03)
inner_2.add_nuclide(‘Pu241’, 2.3422E-04)
inner_2.add_nuclide(‘Pu242’, 3.3784E-04)
inner_2.add_nuclide(‘Am241’, 3.6180E-05)
inner_2.add_nuclide(‘Am242’, 1.1575E-08)
inner_2.add_nuclide(‘Am242’, 7.7557E-07)
inner_2.add_nuclide(‘Am243’, 1.6267E-05)
inner_2.add_nuclide(‘Cm242’, 1.1575E-08)
inner_2.add_nuclide(‘Cm243’, 1.8567E-06)
inner_2.add_nuclide(‘Cm244’, 5.3407E-08)
inner_2.add_nuclide(‘Cm243’, 1.4363E-06)
inner_2.add_nuclide(‘Cm246’, 4.4683E-08)
inner_2.add_element(‘Mo’, 1.3400E-03)
inner_2.set_density(‘sum’)
#40.22 cm - 60.33 cm from bottom active core
inner_3 = openmc.Material(name = ‘40.22 cm - 60.33 cm from bottom active core’)
inner_3.volume = np.pi*(0.34700.3954)20.1146930
inner_3.temperature = fuel_temperature
inner_3.add_nuclide(‘C12’, 2.6385E-02)
inner_3.add_nuclide(‘U234’, 1.1555E-06)
inner_3.add_nuclide(‘U235’, 4.3713E-05)
inner_3.add_nuclide(‘U236’, 3.3299E-06)
inner_3.add_nuclide(‘U238’, 2.2084E-02)
inner_3.add_nuclide(‘Np237’, 3.7208E-06)
inner_3.add_nuclide(‘Np239’, 4.7916E-06)
inner_3.add_nuclide(‘Pu238’, 1.0175E-04)
inner_3.add_nuclide(‘Pu239’, 1.8745E-03)
inner_3.add_nuclide(‘Pu240’, 1.0165E-03)
inner_3.add_nuclide(‘Pu241’, 2.3224E-04)
inner_3.add_nuclide(‘Pu242’, 3.3674E-04)
inner_3.add_nuclide(‘Am241’, 3.5616E-05)
inner_3.add_nuclide(‘Am242’, 1.2270E-08)
inner_3.add_nuclide(‘Am242’, 8.1635E-07)
inner_3.add_nuclide(‘Am243’, 1.7351E-05)
inner_3.add_nuclide(‘Cm242’, 1.2270E-08)
inner_3.add_nuclide(‘Cm243’, 1.9718E-06)
inner_3.add_nuclide(‘Cm244’, 6.0550E-08)
inner_3.add_nuclide(‘Cm245’, 1.6492E-06)
inner_3.add_nuclide(‘Cm246’, 5.5003E-08)
inner_3.add_element(‘Mo’, 1.3400E-03)
inner_3.set_density(‘sum’)
#60.33 cm - 80.44 cm from bottom active core
inner_4 = openmc.Material(name = ‘60.33 - 80.44 cm from bottom active core’)
inner_4.volume = np.pi*(0.34700.3954)20.1146930
inner_4.temperature = fuel_temperature
inner_4.add_nuclide(‘C12’, 2.6385E-02)
inner_4.add_nuclide(‘U234’, 1.1735E-06)
inner_4.add_nuclide(‘U235’, 4.5174E-05)
inner_4.add_nuclide(‘U236’, 3.0230E-06)
inner_4.add_nuclide(‘U238’, 2.2177E-02)
inner_4.add_nuclide(‘Np237’, 3.3324E-06)
inner_4.add_nuclide(‘Np239’, 4.2607E-06)
inner_4.add_nuclide(‘Pu238’, 1.0383E-04)
inner_4.add_nuclide(‘Pu239’, 1.8508E-03)
inner_4.add_nuclide(‘Pu240’, 1.0155E-03)
inner_4.add_nuclide(‘Pu241’, 2.3519E-04)
inner_4.add_nuclide(‘Pu242’, 3.3844E-04)
inner_4.add_nuclide(‘Am241’, 3.6509E-05)
inner_4.add_nuclide(‘Am242’, 1.1142E-08)
inner_4.add_nuclide(‘Am242’, 7.4933E-07)
inner_4.add_nuclide(‘Am243’, 1.5594E-05)
inner_4.add_nuclide(‘Cm242’, 1.1142E-08)
inner_4.add_nuclide(‘Cm243’, 1.7856E-06)
inner_4.add_nuclide(‘Cm244’, 4.9215E-08)
inner_4.add_nuclide(‘Cm245’, 1.3185E-06)
inner_4.add_nuclide(‘Cm246’, 3.9501E-08)
inner_4.add_element(‘Mo’, 1.3400E-03)
inner_4.set_density(‘sum’)
#80.44 - 100.55 cm from bottom active core
inner_5 = openmc.Material(name = ‘80.44 - 100.55 cm from bottom active core’)
inner_5.volume = np.pi*(0.34700.3954)20.1146930
inner_5.temperature = fuel_temperature
inner_5.add_nuclide(‘C12’, 2.6385E-02)
inner_5.add_nuclide(‘U234’, 1.2112E-06)
inner_5.add_nuclide(‘U235’, 4.8320E-05)
inner_5.add_nuclide(‘U236’, 2.4244E-06)
inner_5.add_nuclide(‘U238’, 2.2382E-02)
inner_5.add_nuclide(‘Np237’, 2.2640E-06)
inner_5.add_nuclide(‘Np239’, 3.1154E-06)
inner_5.add_nuclide(‘Pu238’, 1.0844E-04)
inner_5.add_nuclide(‘Pu239’, 1.7942E-03)
inner_5.add_nuclide(‘Pu240’, 1.0185E-03)
inner_5.add_nuclide(‘Pu241’, 2.4351E-04)
inner_5.add_nuclide(‘Pu242’, 3.4252E-04)
inner_5.add_nuclide(‘Am241’, 3.8462E-05)
inner_5.add_nuclide(‘Am242’, 8.7856E-09)
inner_5.add_nuclide(‘Am242’, 6.0460E-07)
inner_5.add_nuclide(‘Am243’, 1.2110E-05)
inner_5.add_nuclide(‘Cm242’, 8.7856E-09)
inner_5.add_nuclide(‘Cm243’, 1.3992E-06)
inner_5.add_nuclide(‘Cm244’, 3.0944E-08)
inner_5.add_nuclide(‘Cm245’, 7.7138E-07)
inner_5.add_nuclide(‘Cm246’,1.7550E-08)
inner_5.add_element(‘Mo’, 1.3400E-03)
inner_5.set_density(‘sum’)
#0 - 20.11 cm from bottom active core
outer_1 = openmc.Material(name = ‘outer 0 - 20.11 cm from bottom active core’)
outer_1.volume = np.pi*(0.34700.3954)20.1146930
outer_1.temperature = fuel_temperature
outer_1.add_nuclide(‘C12’, 2.6385E-02)
outer_1.add_nuclide(‘U234’, 2.2313E-06)
outer_1.add_nuclide(‘U235’, 3.6556E-05)
outer_1.add_nuclide(‘U236’, 3.8682E-06)
outer_1.add_nuclide(‘U238’, 1.9613E-02)
outer_1.add_nuclide(‘Np237’, 4.0477E-06)
outer_1.add_nuclide(‘Np239’, 3.0567E-06)
outer_1.add_nuclide(‘Pu238’, 1.2031E-04)
outer_1.add_nuclide(‘Pu239’, 2.1078E-03)
outer_1.add_nuclide(‘Pu240’, 1.2520E-03)
outer_1.add_nuclide(‘Pu241’, 2.7216E-04)
outer_1.add_nuclide(‘Pu242’, 4.0896E-04)
outer_1.add_nuclide(‘Am241’, 5.2029E-05)
outer_1.add_nuclide(‘Am242’, 1.3235E-08)
outer_1.add_nuclide(‘Am242’, 1.3674E-06)
outer_1.add_nuclide(‘Am243’, 2.6522E-05)
outer_1.add_nuclide(‘Cm242’, 2.3638E-06)
outer_1.add_nuclide(‘Cm243’, 1.0326E-07)
outer_1.add_nuclide(‘Cm244’, 3.2259E-06)
outer_1.add_nuclide(‘Cm245’, 1.3656E-07)
outer_1.add_nuclide(‘Cm246’, 3.0878E-09)
outer_1.add_element(‘Mo’, 1.8258E-03)
outer_1.set_density(‘sum’)
#20.11 - 40.22 cm from bottom active core
outer_2 = openmc.Material(name = ‘outer 20.11 - 40.22 cm from bottom active core’)
outer_2.volume = np.pi*(0.34700.3954)20.1146930
outer_2.temperature = fuel_temperature
outer_2.add_nuclide(‘C12’, 2.6385E-02)
outer_2.add_nuclide(‘U234’, 2.1245E-06)
outer_2.add_nuclide(‘U235’, 3.3005E-05)
outer_2.add_nuclide(‘U236’, 4.4412E-06)
outer_2.add_nuclide(‘U238’, 1.9319E-02)
outer_2.add_nuclide(‘Np237’, 5.6298E-06)
outer_2.add_nuclide(‘Np239’, 3.9393E-06)
outer_2.add_nuclide(‘Pu238’, 1.1288E-04)
outer_2.add_nuclide(‘Pu239’, 2.1362E-03)
outer_2.add_nuclide(‘Pu240’, 1.2386E-03)
outer_2.add_nuclide(‘Pu241’, 2.5863E-04)
outer_2.add_nuclide(‘Pu242’, 4.0010E-04)
outer_2.add_nuclide(‘Am241’, 4.8117E-05)
outer_2.add_nuclide(‘Am242’, 1.5409E-08)
outer_2.add_nuclide(‘Am242’, 1.5394E-06)
outer_2.add_nuclide(‘Am243’, 3.1778E-05)
outer_2.add_nuclide(‘Cm242’, 2.7689E-06)
outer_2.add_nuclide(‘Cm243’, 1.3838E-07)
outer_2.add_nuclide(‘Cm244’, 4.8143E-06)
outer_2.add_nuclide(‘Cm245’, 2.5044E-07)
outer_2.add_nuclide(‘Cm246’, 7.2444E-09)
outer_2.add_element(‘Mo’, 1.8258E-03)
outer_2.set_density(‘sum’)
#40.22 - 60.33 cm from bottom active core
outer_3 = openmc.Material(name = ‘outer 40.22 - 60.33 cm from bottom active core’)
outer_3.volume = np.pi*(0.34700.3954)20.1146930
outer_3.temperature = fuel_temperature
outer_3.add_nuclide(‘C12’, 2.6385E-02)
outer_3.add_nuclide(‘U234’, 2.0828E-06)
outer_3.add_nuclide(‘U235’, 3.1658E-05)
outer_3.add_nuclide(‘U236’, 4.7031E-06)
outer_3.add_nuclide(‘U238’, 1.9211E-02)
outer_3.add_nuclide(‘Np237’, 6.0470E-06)
outer_3.add_nuclide(‘Np239’, 4.2728E-06)
outer_3.add_nuclide(‘Pu238’, 1.1019E-04)
outer_3.add_nuclide(‘Pu239’, 2.1436E-03)
outer_3.add_nuclide(‘Pu240’, 1.2385E-03)
outer_3.add_nuclide(‘Pu241’, 2.5524E-04)
outer_3.add_nuclide(‘Pu242’, 3.9726E-04)
outer_3.add_nuclide(‘Am241’, 4.6701E-05)
outer_3.add_nuclide(‘Am242’, 1.6296E-08)
outer_3.add_nuclide(‘Am242’, 1.6065E-06)
outer_3.add_nuclide(‘Am243’, 3.4010E-05)
outer_3.add_nuclide(‘Cm242’, 2.9338E-06)
outer_3.add_nuclide(‘Cm243’, 1.5726E-07)
outer_3.add_nuclide(‘Cm244’, 5.5901E-06)
outer_3.add_nuclide(‘Cm245’, 3.1218E-07)
outer_3.add_nuclide(‘Cm246’, 9.7807E-09)
outer_3.add_element(‘Mo’, 1.8258E-03)
outer_3.set_density(‘sum’)
#60.33 - 80.44 cm from bottom active core
outer_4 = openmc.Material(name = ‘outer 60.33 - 80.44 cm from bottom active core’)
outer_4.volume = np.pi*(0.34700.3954)20.1146930
outer_4.temperature = fuel_temperature
outer_4.add_nuclide(‘C12’, 2.6385E-02)
outer_4.add_nuclide(‘U234’, 2.1405E-06)
outer_4.add_nuclide(‘U235’, 3.3521E-05)
outer_4.add_nuclide(‘U236’, 4.3381E-06)
outer_4.add_nuclide(‘U238’, 1.9359E-02)
outer_4.add_nuclide(‘Np237’, 5.4769E-06)
outer_4.add_nuclide(‘Np239’, 3.8143E-06)
outer_4.add_nuclide(‘Pu238’, 1.1390E-04)
outer_4.add_nuclide(‘Pu239’, 2.1335E-03)
outer_4.add_nuclide(‘Pu240’, 1.2384E-03)
outer_4.add_nuclide(‘Pu241’, 2.5984E-04)
outer_4.add_nuclide(‘Pu242’, 4.0115E-04)
outer_4.add_nuclide(‘Am241’, 4.8656E-05)
outer_4.add_nuclide(‘Am242’, 1.5069E-08)
outer_4.add_nuclide(‘Am242’, 1.5130E-06)
outer_4.add_nuclide(‘Am243’, 3.0910E-05)
outer_4.add_nuclide(‘Cm242’, 2.7054E-06)
outer_4.add_nuclide(‘Cm243’, 1.3092E-07)
outer_4.add_nuclide(‘Cm244’, 4.5164E-06)
outer_4.add_nuclide(‘Cm245’, 2.2707E-07)
outer_4.add_nuclide(‘Cm246’, 6.3131E-09)
outer_4.add_element(‘Mo’, 1.8258E-03)
outer_4.set_density(‘sum’)
#80.44 - 100.55 cm from bottom active core
outer_5= openmc.Material(name = ‘outer 80.44 - 100.55 cm from bottom active core’)
outer_5.volume = np.pi*(0.34700.3954)20.1146930
outer_5.temperature = fuel_temperature
outer_5.add_nuclide(‘C12’, 2.6385E-02)
outer_5.add_nuclide(‘U234’, 2.2687E-06)
outer_5.add_nuclide(‘U235’, 3.7812E-05)
outer_5.add_nuclide(‘U236’, 3.5850E-06)
outer_5.add_nuclide(‘U238’, 1.9697E-02)
outer_5.add_nuclide(‘Np237’, 3.7682E-06)
outer_5.add_nuclide(‘Np239’, 2.7851E-06)
outer_5.add_nuclide(‘Pu238’, 1.2263E-04)
outer_5.add_nuclide(‘Pu239’, 2.1019E-03)
outer_5.add_nuclide(‘Pu240’, 1.2492E-03)
outer_5.add_nuclide(‘Pu241’, 2.7437E-04)
outer_5.add_nuclide(‘Pu242’, 4.1086E-04)
outer_5.add_nuclide(‘Am241’, 5.3279E-05)
outer_5.add_nuclide(‘Am242’, 1.2281E-08)
outer_5.add_nuclide(‘Am242’, 1.2836E-06)
outer_5.add_nuclide(‘Am243’, 2.4305E-05)
outer_5.add_nuclide(‘Cm242’, 2.1894E-06)
outer_5.add_nuclide(‘Cm243’, 8.5277E-08)
outer_5.add_nuclide(‘Cm244’, 2.6306E-06)
outer_5.add_nuclide(‘Cm245’, 9.9197E-08)
outer_5.add_nuclide(‘Cm246’, 1.9680E-09)
outer_5.add_element(‘Mo’, 1.8258E-03)
outer_5.set_density(‘sum’)
#spesifikasi material cladding, reflector, pendingin, absorber, dan lain-lain
em10_steel = openmc.Material(name = ‘Lower and upper structure’)
em10_steel.temperature = structure_temperature
em10_steel.add_element(‘C’, 3.8254E-04)
em10_steel.add_element(‘Si’, 4.9089E-04)
em10_steel.add_element(‘Ti’, 1.9203E-05)
em10_steel.add_element(‘Cr’, 7.5122E-03)
em10_steel.add_element(‘Fe’, 7.3230E-02)
em10_steel.add_element(‘Ni’, 3.9162E-04)
em10_steel.add_element(‘Mo’, 4.7925E-04)
em10_steel.add_nuclide(‘Mn55’, 4.1817E-04)
em10_steel.set_density(‘sum’)
sodium = openmc.Material(name=‘pendingin’)
sodium.temperature = structure_temperature
sodium.add_element(‘Na’, 2.1924E-02)
sodium.set_density(‘sum’)
clad = openmc.Material(name=‘Cladding atau kelongsong’)
clad.temperature = structure_temperature
clad.add_element(‘C’, 3.5740E-04)
clad.add_element(‘O’, 3.9924E-04)
clad.add_element(‘Ti’, 5.3824E-04)
clad.add_element(‘Cr’, 1.7753E-02)
clad.add_element(‘Fe’, 5.3872E-02)
clad.add_element(‘Ni’, 3.6588E-04)
clad.add_nuclide(‘Mn55’, 2.3441E-04)
clad.add_element(‘P’, 2.7718E-05)
clad.add_element(‘Al’, 9.1482E-03)
clad.add_element(‘Co’, 2.1852E-04)
clad.add_element(‘Cu’, 1.0135E-04)
clad.add_element(‘Y’, 2.6616E-04)
nat_b4c = openmc.Material(name=‘natural boron karbida’)
nat_b4c.temperature = structure_temperature
nat_b4c.add_element(‘C’,0.00027)
nat_b4c.add_nuclide(‘B10’, 0.00981)
nat_b4c.add_nuclide(‘B11’, 0.00991)
nat_b4c.set_density(‘sum’)
‘’’
en_b4c = openmc.Material(name=‘boron karbida diperkaya’)
en_b4c.temperature = structure_temperature
en_b4c.add_element(‘C’,0.020632)
en_b4c.add_nuclide(‘B10’, 0.053642)
en_b4c.add_nuclide(‘B11’, 0.028884)
en_b4c.set_density(‘sum’)
‘’’
helium = openmc.Material(name=‘helium untuk gas plenum’)
helium.temperature = structure_temperature
helium.set_density(‘g/cm3’, 0.001598)
helium.add_element(‘He’, 2.4044E-4)
primary_control = openmc.Material(name=‘primary control rod’)
primary_control.temperature = structure_temperature
primary_control.add_element(‘C’, 2.70E-02)
primary_control.add_nuclide(‘B10’, 2.3E-02)
primary_control.add_nuclide(‘B11’, 8.4E-02)
primary_control.set_density(‘sum’)
sec_control = openmc.Material(name=‘secondary control rod’)
sec_control.temperature = structure_temperature
sec_control.add_element(‘C’, 2.70E-02)
sec_control.add_nuclide(‘B10’, 9.8E-02)
sec_control.add_nuclide(‘B11’, 9.9E-02)
sec_control.set_density(‘sum’)
#Export Material XML
materials_file = openmc.Materials([inner_1, inner_2, inner_3, inner_4, inner_5, outer_1, outer_2, outer_3, outer_4, outer_5, em10_steel, sodium, clad, helium])
materials_file.export_to_xml()
#Fill Empty Space
Membuat isi dari space kosong dengan sodium
sodium_mod_cell = openmc.Cell(cell_id=1, fill=sodium)
sodium_mod_u = openmc.Universe(universe_id=1, cells=(sodium_mod_cell,))
#surface and region
fuel_radius_inner = openmc.ZCylinder(r=0.3319)
clad_inner_radius_inner = openmc.ZCylinder(r=0.3470)
clad_outer_radius_inner = openmc.ZCylinder(r=0.3954)
sleeve_min_z = openmc.ZPlane(z0=-50, boundary_type = ‘transmission’)
sleeve_max_z = openmc.ZPlane(z0=+50, boundary_type = ‘transmission’)
axial_reflector_max_z = openmc.ZPlane(z0= +131, boundary_type = ‘vacuum’)
axial_reflector_min_z = openmc.ZPlane(z0= -80, boundary_type = ‘transmission’)
gas_plenum_max_z = openmc.ZPlane(z0 = +61, boundary_type = ‘transmission’)
gas_plenum_min_z = openmc.ZPlane (z0 = -171, boundary_type = ‘vacuum’)
fuel_region_inner = -fuel_radius_inner & -sleeve_max_z & +sleeve_min_z
gap_region_inner = +fuel_radius_inner & -clad_inner_radius_inner & -sleeve_max_z & +sleeve_min_z
clad_region_inner = +clad_inner_radius_inner & -clad_outer_radius_inner & -axial_reflector_max_z & +gas_plenum_min_z
axial_reflector_region_top_inner = -clad_inner_radius_inner & -axial_reflector_max_z & +gas_plenum_max_z
axial_reflector_region_bottom_inner = -clad_inner_radius_inner & +axial_reflector_min_z & -sleeve_min_z
gas_plenum_region_top_inner = -clad_inner_radius_inner & +sleeve_max_z & -gas_plenum_max_z
gas_plenum_region_bottom_inner = -clad_inner_radius_inner & +gas_plenum_min_z & -axial_reflector_min_z
surrounding_region_inner = +clad_outer_radius_inner & -axial_reflector_max_z & +gas_plenum_min_z
fuell_cell = openmc.Cell(fill=inner_1, region=fuel_region_inner)
gap_cell_inner = openmc.Cell(fill=helium, region=gap_region_inner)
clad_cell_inner = openmc.Cell(fill=clad, region=clad_region_inner)
axial_reflector_top_cell_inner = openmc.Cell(fill=sodium, region=axial_reflector_region_top_inner)
axial_reflector_bottom_cell_inner = openmc.Cell(fill=sodium, region=axial_reflector_region_bottom_inner)
gas_plenum_cell_top_inner = openmc.Cell(fill=helium, region=gas_plenum_region_top_inner)
gas_plenum_cell_bottom_inner = openmc.Cell(fill=helium, region=gas_plenum_region_bottom_inner)
surrounding_cell_inner = openmc.Cell(fill=sodium, region=surrounding_region_inner)
inner_universe = openmc.Universe(cells=(fuell_cell, gap_cell_inner, clad_cell_inner, axial_reflector_top_cell_inner, axial_reflector_bottom_cell_inner, gas_plenum_cell_top_inner, gas_plenum_cell_bottom_inner, surrounding_cell_inner))
#hexagonal region
sodium_fill_cell = openmc.Cell(fill=sodium)
sodium_universe = openmc.Universe(cells=(sodium_fill_cell,))
sodium_prism = openmc.hexagonal_prism(edge_length=0.99, orientation=‘y’)
sodium_prism_region = sodium_prism & +gas_plenum_min_z & -axial_reflector_max_z
sodium_prism_out = ~sodium_prism & +gas_plenum_min_z & -axial_reflector_max_z
sodium_prism_cell = openmc.Cell(fill=sodium, region=sodium_prism_region)
sodium_prism_cell_out = openmc.Cell(fill=sodium, region=sodium_prism_out)
sodium_prism_universe = openmc.Universe(cells=(sodium_prism_cell, sodium_prism_cell_out))
#inner assemblies
inner_lattice = openmc.HexLattice()
inner_lattice.center = (0., 0.)
inner_lattice.pitch = (21.08/17,)
inner_lattice.orientation = ‘y’
inner_lattice.outer = sodium_universe
inner_ring_one = [inner_universe] * 48
inner_ring_two = [inner_universe] * 42
inner_ring_three = [inner_universe] * 36
inner_ring_four = [inner_universe] * 30
inner_ring_five = [inner_universe] * 24
inner_ring_six = [inner_universe] * 18
inner_ring_seven = [inner_universe] * 12
inner_ring_eight = [inner_universe] * 6
inner_ring_nine = [inner_universe]
inner_lattice.universes = [inner_ring_one, inner_ring_two, inner_ring_three, inner_ring_four, inner_ring_five, inner_ring_six, inner_ring_seven, inner_ring_eight, inner_ring_nine]
inner_prism = openmc.model.hexagonal_prism(edge_length=9.3803, orientation=‘x’)
wrapper_tube_inner = openmc.model.hexagonal_prism(edge_length=9.3803, orientation=‘x’)
wrapper_tube_region_inner = (~inner_prism & wrapper_tube_inner & +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_inner = (~wrapper_tube_inner & -axial_reflector_max_z & +gas_plenum_min_z)
main_in_assembly_cell = openmc.Cell(fill=inner_lattice, region=inner_prism & -axial_reflector_max_z & +gas_plenum_min_z)
wrapper_tube_inner_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_inner)
wrapper_tube_inner_outside_cell = openmc.Cell(fill=sodium, region=outside_wrapper_tube_inner)
inner_assembly_universe = openmc.Universe(cells=[main_in_assembly_cell, wrapper_tube_inner_cell, wrapper_tube_inner_outside_cell])
#####outer pincell
fuel_radius_outer = openmc.ZCylinder(r=0.3319)
clad_inner_radius_outer = openmc.ZCylinder(r=0.3470)
clad_outer_radius_outer = openmc.ZCylinder(r=0.3954)
fuel_region_outer = -fuel_radius_outer & -sleeve_max_z & +sleeve_min_z
gap_region_outer = +fuel_radius_outer & -clad_inner_radius_outer & -sleeve_max_z & +sleeve_min_z
clad_region_outer = +clad_inner_radius_outer & -clad_outer_radius_outer & -axial_reflector_max_z & +gas_plenum_min_z
axial_reflector_region_top_outer = -clad_inner_radius_outer & -axial_reflector_max_z & +gas_plenum_max_z
axial_reflector_region_bottom_outer = -clad_inner_radius_outer & +axial_reflector_min_z & -sleeve_min_z
gas_plenum_region_top_outer = -clad_inner_radius_outer & +sleeve_max_z & -gas_plenum_max_z
gas_plenum_region_bottom_outer = -clad_inner_radius_outer & +gas_plenum_min_z & -axial_reflector_min_z
surrounding_region_outer = +clad_outer_radius_outer & -axial_reflector_max_z & +gas_plenum_min_z
outer_fuel_cell = openmc.Cell(fill=outer_1, region=fuel_region_outer)
gap_cell_outer = openmc.Cell(fill=helium, region=gap_region_outer)
clad_cell_outer = openmc.Cell(fill=clad, region=clad_region_outer)
axial_reflector_top_cell_outer = openmc.Cell(fill=sodium, region=axial_reflector_region_top_outer)
axial_reflector_bottom_cell_outer = openmc.Cell(fill=sodium, region=axial_reflector_region_bottom_outer)
gas_plenum_cell_top_outer = openmc.Cell(fill = helium, region=gas_plenum_region_top_outer)
gas_plenum_cell_bottom_outer = openmc.Cell(fill = helium, region=gas_plenum_region_bottom_outer)
surrounding_cell_outer = openmc.Cell(fill=sodium, region=surrounding_region_outer)
outer_universe = openmc.Universe(cells=(outer_fuel_cell, gap_cell_outer, clad_cell_outer, axial_reflector_top_cell_outer, axial_reflector_bottom_cell_outer, gas_plenum_cell_top_outer, gas_plenum_cell_bottom_outer, surrounding_cell_outer))
#outer assemblies
outer_lattice = openmc.HexLattice()
outer_lattice.center = (0., 0.)
outer_lattice.pitch = (21.08/17,)
outer_lattice.orientation = ‘x’
outer_lattice.outer = sodium_universe
out_ring_zero = [outer_universe] * 54
outer_ring_one = [outer_universe] * 48
outer_ring_two = [outer_universe] * 42
outer_ring_three = [outer_universe] * 36
outer_ring_four = [outer_universe] * 30
outer_ring_five = [outer_universe] * 24
outer_ring_six = [outer_universe] * 18
outer_ring_seven = [outer_universe] * 12
outer_ring_eight = [outer_universe] * 6
outer_ring_nine = [outer_universe]
outer_lattice.universes = [out_ring_zero, outer_ring_one, outer_ring_two, outer_ring_three, outer_ring_four, outer_ring_five, outer_ring_six, outer_ring_seven, outer_ring_eight, outer_ring_nine]
outer_prism = openmc.model.hexagonal_prism(edge_length=9.3803, orientation=‘x’)
wrapper_tube_outer = openmc.model.hexagonal_prism(edge_length=9.3803, orientation=‘x’)
wrapper_tube_region_outer = (~outer_prism & wrapper_tube_outer & +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_outer = (~wrapper_tube_outer & -axial_reflector_max_z & +gas_plenum_min_z)
main_out_assembly_cell = openmc.Cell(fill=outer_lattice, region=outer_prism & -axial_reflector_max_z & +gas_plenum_min_z)
wrapper_tube_outer_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_outer)
wrapper_tube_outer_outside_cell = openmc.Cell(fill=sodium, region=outside_wrapper_tube_outer)
outer_assembly_universe = openmc.Universe(cells=[main_out_assembly_cell, wrapper_tube_outer_cell, wrapper_tube_outer_outside_cell])
#reflector region
reflector_radius = openmc.ZCylinder(r=0.3319)
reflector_cladding_radius = openmc.ZCylinder(r=0.7424)
reflector_region = -reflector_radius & +gas_plenum_min_z & -axial_reflector_max_z
reflector_clad_region = (+reflector_radius & -reflector_cladding_radius & +gas_plenum_min_z & -axial_reflector_max_z)
reflector_cell = openmc.Cell(fill=sodium, region=reflector_region)
reflector_clad_cell = openmc.Cell(fill=em10_steel, region=reflector_clad_region)
reflector_out_cell = openmc.Cell(fill=sodium, region = +reflector_cladding_radius & +gas_plenum_min_z & -axial_reflector_max_z)
reflector_universe = openmc.Universe(cells = (reflector_cell, reflector_clad_cell, reflector_out_cell))
‘’’
#control rod
control_rod_radius = openmc.ZCylinder(r=0.915)
control_clad_inner_radius = openmc.ZCylinder(r=1.0415)
control_clad_outer_radius = openmc.ZCylinder(r=1.1412)
control_rod_region = (-control_rod_radius & +sleeve_max_z & -axial_reflector_max_z)
control_gap_region = (+control_rod_radius & -control_clad_inner_radius & +sleeve_max_z & -axial_reflector_max_z)
control_clad_region = (+control_clad_inner_radius & -control_clad_outer_radius & +gas_plenum_min_z & -axial_reflector_max_z)
follower_region = (-control_clad_outer_radius & +gas_plenum_min_z & -sleeve_max_z)
control_outside_region = (+control_clad_outer_radius & +gas_plenum_min_z & -axial_reflector_max_z)
control_rod_cell = openmc.Cell(fill=sodium, region=control_rod_region)
control_gap_cell = openmc.Cell(fill=helium, region=control_gap_region)
control_clad_cell = openmc.Cell(fill=em10_steel, region=control_clad_region)
follower_cell = openmc.Cell(fill=sodium, region=follower_region)
control_surroundings_cell = openmc.Cell(fill=sodium, region=control_outside_region)
control_universe = openmc.Universe(cells=(control_rod_cell, control_gap_cell, control_clad_cell, follower_cell, control_surroundings_cell))
‘’’
#reflector assemblies
reflector_lattice = openmc.HexLattice()
reflector_lattice.center = (0., 0.)
reflector_lattice.pitch = (21.08/17,)
reflector_lattice.orientation = ‘x’
reflector_lattice.outer = sodium_universe
reflector_inner_ring_one = [reflector_universe] * 12
reflector_inner_ring_two = [reflector_universe] * 6
reflector_inner_ring_three = [reflector_universe]
reflector_lattice.universes = [reflector_inner_ring_one, reflector_inner_ring_two, reflector_inner_ring_three]
reflector_prism = openmc.model.hexagonal_prism(edge_length=11.65092843, orientation=‘x’)
reflector_wrapper =openmc.model.hexagonal_prism(edge_length=11.91073605, orientation=‘x’)
wrapper_tube_region_ref = (~reflector_prism & reflector_wrapper & +gas_plenum_min_z & -axial_reflector_max_z)
outside_wrapper_tube_ref = (~reflector_wrapper & -axial_reflector_max_z & +gas_plenum_min_z)
wrapper_tube_ref_cell = openmc.Cell(fill=em10_steel, region=wrapper_tube_region_ref)
wrapper_tube_outside_ref_cell = openmc.Cell(fill=sodium, region=outside_wrapper_tube_ref)
main_reflector_assembly_cell = openmc.Cell(fill=reflector_lattice, region=reflector_prism & -axial_reflector_max_z & +gas_plenum_min_z)
reflector_assembly_universe = openmc.Universe(cells=[main_reflector_assembly_cell, wrapper_tube_ref_cell, wrapper_tube_outside_ref_cell])
#core lattice
main_lattice = openmc.HexLattice()
main_lattice.center = (0., 0.)
main_lattice.pitch = (21.08,)
main_lattice.orientation = ‘y’ #barudiubah
main_lattice.outer = sodium_universe
#fully up
reflector_ring_one = [reflector_assembly_universe] * 96
reflector_ring_two = [reflector_assembly_universe] * 90
reflector_ring_three = [reflector_assembly_universe] * 84
reflector_ring_four = ([reflector_assembly_universe] * 5 + [outer_assembly_universe] * 4 + [reflector_assembly_universe] * 4) * 6
reflector_ring_five = ([reflector_assembly_universe] + [outer_assembly_universe] * 11) * 6
out_ring_one = [outer_assembly_universe] * 66
out_ring_two = ([outer_assembly_universe] * 5 + [inner_assembly_universe] + ([outer_assembly_universe] * 9 + [inner_assembly_universe]) * 5 + [outer_assembly_universe] * 4)
out_ring_three = ([outer_assembly_universe] * 2 + ([inner_assembly_universe] + [inner_assembly_universe] * 4 + [inner_assembly_universe] + [outer_assembly_universe] * 3) * 5 + [inner_assembly_universe] + [inner_assembly_universe] * 4 + [inner_assembly_universe] + [outer_assembly_universe])
in_ring_one = [inner_assembly_universe] * 48
in_ring_two = [inner_assembly_universe] * 42
in_ring_three = ([inner_assembly_universe]*3 + [sodium_prism_universe]) * 9
in_ring_four = [inner_assembly_universe] * 30
in_ring_five = [inner_assembly_universe] * 24
in_ring_six = ([inner_assembly_universe] + [inner_assembly_universe] * 2) * 6
in_ring_seven = [inner_assembly_universe] * 12
in_ring_eight = [inner_assembly_universe] * 6
in_ring_nine = [reflector_assembly_universe] * 1
main_lattice.universes = [reflector_ring_one, reflector_ring_two, reflector_ring_three, reflector_ring_four, reflector_ring_five, out_ring_one,out_ring_two,out_ring_three,in_ring_one,in_ring_two, in_ring_three,in_ring_four,in_ring_five,in_ring_six, in_ring_seven,in_ring_eight,in_ring_nine]
core_surface = openmc.model.hexagonal_prism(edge_length=347.82, boundary_type=‘vacuum’)
core = openmc.Cell(fill=main_lattice, region = core_surface & -axial_reflector_max_z & +gas_plenum_min_z)
outside_core = openmc.Cell(fill=sodium, region = ~core_surface & -axial_reflector_max_z & +gas_plenum_min_z)
inside_core = openmc.Cell(fill=sodium, region = ~core_surface & -axial_reflector_max_z & +gas_plenum_min_z)
core_universe = openmc.Universe(universe_id=16, cells=[core, ])
#geometry xml
geom = openmc.Geometry(core_universe)
geom.export_to_xml()
####### PLOT XML ####################
Plot dengan orientasi xy
plot_xy = openmc.Plot()
plot_xy.origin = (0, 0, 0)
plot_xy.width = (800, 800)
plot_xy.pixels = (3000, 3000)
plot_xy.filename = ‘CORE XY’
plot_xy.color_by = ‘cell’
plot_xy.basis = ‘xy’
plots_xy = openmc.Plots([plot_xy])
plots_xy.export_to_xml()
openmc.plot_geometry()
Plot dengan orientasi xz
plot_xz = openmc.Plot()
plot_xz.origin = (0, 0, 0)
plot_xz.width = (100, 50)
plot_xz.pixels = (1000, 1000)
plot_xz.filename = ‘CORE XZ’
plot_xz.color_by = ‘cell’
plot_xz.basis = ‘xz’
plots_xz = openmc.Plots([plot_xz])
plots_xz.export_to_xml()
openmc.plot_geometry()
Plot dengan orientasi yz
plot_yz = openmc.Plot()
plot_yz.origin = (0, 0, 0)
plot_yz.width = (100, 50)
plot_yz.pixels = (1000, 1000)
plot_yz.filename = ‘CORE YZ’
plot_yz.color_by = ‘cell’
plot_yz.basis = ‘yz’
plots_yz = openmc.Plots([plot_yz])
plots_yz.export_to_xml()
openmc.plot_geometry()
Plot dengan orientasi zy
plot_zy = openmc.Plot()
plot_zy.origin = (0, 0, 0)
plot_zy.width = (800, 800)
plot_zy.pixels = (3000, 3000)
plot_zy.filename = ‘CORE ZY’
plot_zy.color_by = ‘cell’
plots_zy = openmc.Plots([plot_zy])
plots_zy.export_to_xml()
openmc.plot_geometry()
duct_thickness = 0.4525
subassembly_duct = 20.2641
pitch = 20.9889
PEMBUATAN TALLY XML
Define surfaces used to construct regions
H = 301.70
D = pitch*15
zmin, zmax, radius = -H/2, H/2, D/2
#energy filter
energy_filter = openmc.EnergyFilter([3.6e6,10e6])
#univ filter
univ_filter = openmc.UniverseFilter(inner_assembly_universe)
#Spatial Legendre Filter
Legendre_filter = openmc.SpatialLegendreFilter(50, ‘z’, zmin, zmax)
#Zernike Filter
Zernike_filter = openmc.ZernikeFilter(order=50, x=0.0, y=0.0, r=radius)
#Zernike Radial Filter
Zernikeradial_filter = openmc.ZernikeRadialFilter(order=50, x=0.0, y=0.0, r=radius)
Create a fission tally using Legendre spatial filter
fis_tally_Legendre = openmc.Tally(name = ‘fission Legendre’)
fis_tally_Legendre.scores = [‘fission’]
fis_tally_Legendre.filters = [univ_filter, Legendre_filter]
Create a Zernike azimuthal polinomial expansion filter and add to tally
fis_tally_Zernike = openmc.Tally(name = ‘fission Zernike’)
fis_tally_Zernike.scores = [‘fission’]
fis_tally_Zernike.filters = [univ_filter, Zernike_filter]
Create a Zernike radial polinomial expansion filter and add to tally
fis_tally_Zernikeradial = openmc.Tally(name = ‘fission Zernike radial’)
fis_tally_Zernikeradial.scores = [‘fission’]
fis_tally_Zernikeradial.filters = [univ_filter, Zernikeradial_filter]
Create a flux tally using Legendre spatial filter
flux_tally_Legendre = openmc.Tally(name = ‘flux Legendre’)
flux_tally_Legendre.scores = [‘flux’]
flux_tally_Legendre.filters = [univ_filter, Legendre_filter]
Create a Zernike azimuthal polinomial expansion filter and add to tally
flux_tally_Zernike = openmc.Tally(name = ‘flux Zernike’)
flux_tally_Zernike.scores = [‘flux’]
flux_tally_Zernike.filters = [univ_filter, Zernike_filter]
Create a Zernike radial polinomial expansion filter and add to tally
flux_tally_Zernikeradial = openmc.Tally(name = ‘flux Zernike radial’)
flux_tally_Zernikeradial.scores = [‘flux’]
flux_tally_Zernikeradial.filters = [univ_filter, Zernikeradial_filter]
Heating Tally
heating_tally = openmc.Tally(name = ‘heating’)
heating_tally.scores = [‘heating’]
heating_tally.filters = [univ_filter]
#Create tallies
tallies = openmc.Tallies([heating_tally,flux_tally_Legendre,flux_tally_Zernike, flux_tally_Zernikeradial, fis_tally_Legendre, fis_tally_Zernike, fis_tally_Zernikeradial])
tallies.export_to_xml()
####### SETTINGS XML ####################
H = 301.70 ####overal length of subassembly dalam satuan cm
D = pitch*15
source = openmc.Source()
lower_left = [-D/2, -D/2, -H/2]
upper_right = [D/2, D/2, H/2]
source.space = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
settings = openmc.Settings()
settings.temperature= {‘method’ : ‘interpolation’}
settings.source = source
settings.batches = 250#250
settings.inactive = 25 #25
settings.particles = 10000 #Partikelnya bisa diubah-ubah dan partikel ini menentukan lamanya running
settings.export_to_xml()
####### MODEL ####################
model = openmc.model.Model(geometry=geom, settings=settings, tallies=tallies)
model.run()
####### DEPLETION ####################
import openmc.deplete
power = 3.6e9 #ev
chain = openmc.deplete.Chain.from_xml(“./chain_casl_sfr.xml”)
#operator = openmc.deplete.Operator(geom,settings, “./chain_casl_sfr.xml”) #openmc 0.12
operator = openmc.deplete.Operator(model, “./chain_casl_sfr.xml”) #openmc 0.13
max_step = 2operator.heavy_metal/power1e3
timesteps = [28.5,30,30,30,30,30,30,30,30,30,30] #1 tahun
integrator = openmc.deplete.PredictorIntegrator(operator,timesteps,timestep_units =‘d’, power=power)
integrator.integrate()