# Candu depletion problem

**URL:** <https://openmc.discourse.group/t/candu-depletion-problem/3756>\
**Category:** User Support\
**Created:** [December 17, 2023, 9:45pm UTC](https://openmc.discourse.group/t/candu-depletion-problem/3756 "2023-12-17T21:45:38Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![pfproton](https://avatars.discourse-cdn.com/v4/letter/p/d07c76/32.png) [@pfproton](https://openmc.discourse.group/u/pfproton)\
**Post date:** [December 17, 2023, 9:45pm UTC](https://openmc.discourse.group/t/candu-depletion-problem/3756/1 "2023-12-17T21:45:38Z")

</div>

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?

```auto
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()

```