Dear experts, I am trying to make a depletion example for the Candu reactor. I see the k_eff value around 0.4. What adjustments should be made in this Python code to get this value higher?
from math import pi, sin, cos
import numpy as np
import openmc
import json
import numpy as np
from pathlib import Path
import openmc.deplete
# Fuel and clad
temperature_fuel = 960 #Fuel at 960 K with a density of 10.52 g/cc
temperature_structure = 345.15 #Moderator at 345.15 K with a density of 1.086 g/cc
r_fuel = 0.6122
r_clad = 0.6540
pressure_tube_ir = 5.16890
pressure_tube_or = 5.60320
calendria_ir = 6.44780
calendria_or = 6.58750
# Fuel
fuel = openmc.Material(name='fuel', temperature=temperature_fuel)
fuel.add_element('U', 1.0)
fuel.add_element('O', 2.0)
fuel.set_density('g/cm3', 10.0)
fuel.volume = 37 * pi * r_fuel**2
fuel.set_density('g/cm3', 10.52)
# Coolant Water
heavy_water_coolant = openmc.Material(name='heavy water coolant')
heavy_water_coolant.add_nuclide('H2', 0.799449, 'ao' )
heavy_water_coolant.add_nuclide('O16', 0.199768, 'ao' )
heavy_water_coolant.add_nuclide('H1', 0.000783774, 'ao' )
heavy_water_coolant.set_density('g/cm3', 0.806)
# Calandria tube
clad_calandria = openmc.Material(name='Calandria')
clad_calandria.add_element('Zr', 0.160 ) # 16.0%
clad_calandria.add_element('Fe', 0.060 ) # 6.0%
clad_calandria.add_element('Cr', 0.110 ) # 11.0%
clad_calandria.add_element('Ni', 0.9971 ) # 99.71%
clad_calandria.add_element('B', 5.7409e-05 )
clad_calandria.add_element('Li', 2.5259E-04 )
clad_calandria.set_density('g/cm3', 6.4)
# Clad
clad = openmc.Material(name='Clad')
clad.add_element('Zr', 0.160 ) # 16.0%
clad.add_element('Fe', 0.060 ) # 6.0%
clad.add_element('Cr', 0.110 ) # 11.0%
clad.add_element('Ni', 0.9971 ) # 99.71%
clad.add_element('B', 5.7409e-05 )
clad.add_element('Li', 2.5259E-04 )
clad.set_density('g/cm3', 6.3)
# Moder
heavy_water = openmc.Material(name='heavy water')
heavy_water.add_nuclide('H2', 0.798895, 'ao' )
heavy_water.add_nuclide('O16', 0.201016, 'ao' )
heavy_water.add_nuclide('H1', 8.96e-5, 'ao' )
heavy_water.set_density('g/cm3', 1.086)
heavy_water.add_s_alpha_beta('c_D_in_D2O')
# Pressure Tube
tube = openmc.Material(name='PressTube')
tube.add_element('Zr', 0.975 ) # 97.5%
tube.add_element('B', 3.8889E-05 )
tube.add_element('Li', 1.7111E-04 )
tube.set_density('g/cm3', 6.3)
# Radius to center of each ring of fuel pins
ring_radii = np.array([0.0, 1.4885, 2.8755, 4.3305])
# These are the surfaces that will divide each of the rings
radial_surf = [openmc.ZCylinder(R=r) for r in
(ring_radii[:-1] + ring_radii[1:])/2]
water_cells = []
for i in range(ring_radii.size):
# Create annular region
if i == 0:
water_region = -radial_surf[i]
elif i == ring_radii.size - 1:
water_region = +radial_surf[i-1]
else:
water_region = +radial_surf[i-1] & -radial_surf[i]
water_cells.append(openmc.Cell(fill=heavy_water, region=water_region))
bundle_universe = openmc.Universe(cells=water_cells)
surf_fuel = openmc.ZCylinder(R=r_fuel)
surf_clad = openmc.ZCylinder(R=r_clad)
fuel_cell = openmc.Cell(fill=fuel, region=-surf_fuel)
clad_cell = openmc.Cell(fill=clad, region=+surf_fuel & -surf_clad)
cool_cell = openmc.Cell(fill=heavy_water_coolant, region=+surf_clad & -surf_clad)
pin_universe = openmc.Universe(cells=(fuel_cell, clad_cell, cool_cell))
num_pins = [1, 6, 12, 18]
angles = [0, 0, 15, 0]
for i, (r, n, a) in enumerate(zip(ring_radii, num_pins, angles)):
for j in range(n):
# Determine location of center of pin
theta = (a + j/n*360.) * pi/180.
x = r*cos(theta)
y = r*sin(theta)
pin_boundary = openmc.ZCylinder(x0=x, y0=y, R=r_clad)
water_cells[i].region &= +pin_boundary
# That we can identify the pin later when looking at tallies
pin = openmc.Cell(fill=pin_universe, region=-pin_boundary)
pin.translation = (x, y, 0)
pin.id = (i + 1)*100 + j
bundle_universe.add_cell(pin)
pt_inner = openmc.ZCylinder(R=pressure_tube_ir)
pt_outer = openmc.ZCylinder(R=pressure_tube_or)
calendria_inner = openmc.ZCylinder(R=calendria_ir)
calendria_outer = openmc.ZCylinder(R=calendria_or,boundary_type='reflective')
bundle = openmc.Cell(fill=bundle_universe, region=-pt_inner)
pressure_tube = openmc.Cell(fill=tube, region=+pt_inner & -pt_outer)
v1 = openmc.Cell(region=+pt_outer & -calendria_inner)
calendria = openmc.Cell(fill=clad_calandria, region=+calendria_inner & -calendria_outer)
root_universe = openmc.Universe(cells=[bundle, pressure_tube, v1, calendria])
geom = openmc.Geometry(root_universe)
geom.export_to_xml()
mats = openmc.Materials(geom.get_all_materials().values())
mats.export_to_xml()
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 50
settings.inactive = 10
settings.source = openmc.Source(space=openmc.stats.Point())
settings.temperature = {
'default': temperature_structure,
'method': 'interpolation',
}
settings.export_to_xml()
model = openmc.Model(geometry=geom, materials=mats, settings=settings)
chain = openmc.deplete.Chain.from_xml("/mnt/d/openmc/examples/candu/chain_simple.xml")
operator = openmc.deplete.CoupledOperator(model, "/mnt/d/openmc/examples/candu/chain_simple.xml")
power = 174
time_steps = [30] * 6
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power, timestep_units='d')
integrator.integrate()