Candu depletion compare problem

User Support
1

Dear experts, I want to calculate depletion in the candu reactor example below, but when I compared the K_eff result values with serpent, I saw that there were differences.While Serpent K_eff results are around 1, Openmc shows lower values starting from 0.7.How can I solve this problem?

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
# Outer radius of fuel and clad
r_fuel = 0.6122
r_clad = 0.6540
fuel = openmc.Material(name='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

clad = openmc.Material(name='zircaloy')
clad.add_element('Zr', 1.0)
clad.set_density('g/cm3', 6.0)

heavy_water = openmc.Material(name='heavy water')
heavy_water.add_nuclide('H2', 2.0)
heavy_water.add_nuclide('O16', 1.0)
heavy_water.add_s_alpha_beta('c_D_in_D2O')
heavy_water.set_density('g/cm3', 1.1)
# Pressure tube and calendria radii
pressure_tube_ir = 5.16890
pressure_tube_or = 5.60320
calendria_ir = 6.44780
calendria_or = 6.58750

# 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)

fuel_cell = openmc.Cell(fill=fuel, region=-surf_fuel)
clad_cell = openmc.Cell(fill=clad, region=+surf_fuel)

pin_universe = openmc.Universe(cells=(fuel_cell, clad_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
        
        # Create each fuel pin -- note that we explicitly assign an ID so 
        # 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='vacuum')

bundle = openmc.Cell(fill=bundle_universe, region=-pt_inner)
pressure_tube = openmc.Cell(fill=clad, region=+pt_inner & -pt_outer)
v1 = openmc.Cell(region=+pt_outer & -calendria_inner)
calendria = openmc.Cell(fill=clad, region=+calendria_inner & -calendria_outer)

root_universe = openmc.Universe(cells=[bundle, pressure_tube, v1, calendria])
geometry = openmc.Geometry(root_universe)
geometry.export_to_xml()

materials = openmc.Materials([fuel, clad, heavy_water])
materials.export_to_xml()
################################################################################
# SIMULATION SETTINGS

settings = openmc.Settings()
settings.particles = 4_000_000
settings.inactive = 20
settings.batches = 220
settings.statepoint = {'batches': []}
settings.output = {'tallies': False}
settings.export_to_xml()

################################################################################
# create a Model object
model = openmc.Model(geometry=geometry, materials=materials, settings=settings)
# DEPLETION

cumulative_days = np.linspace(0.0, 360.0, 19)
timesteps = np.diff(cumulative_days)
power = 63143.  # 271 * Average linear power from ANL-AFCI-02, Table II.1-3

data_path = Path('/home/user1/openmc/examples/depletion-comparison/data/depletion')
file_path = data_path / 'serpent_fissq.json'

with file_path.open() as fh:
    serpent_q = json.load(fh)


# Create depletion operator
chain_file = data_path / 'chain_simple.xml'
op = openmc.deplete.CoupledOperator(model, chain_file,
    fission_q=serpent_q,
    fission_yield_mode="average"
)

# Execute depletion using integrator
integrator = openmc.deplete.PredictorIntegrator(
    op, timesteps, power, timestep_units='d')
integrator.integrate()
2

Hi. I’ve had a very quick look at the model problem you’ve provided. I don’t see an obvious reason for the OpenMC to SERPENT differences you’ve noted. CNL has modeled the CANDU problem extensively using MCNP, SERPENT, KENO, and OpenMC. We have never observed differences of the magnitude you’re reporting, which suggests an error in the OpenMC model.

Are the material temperatures and densities used in your example the same as those used in the SERPENT model? If not, I would start by making those consistent. They are, by the way, inconsistent with nominal CANDU operating conditions. I would recommend the following conditions:
Fuel at 960 K with a density of 10.52 g/cc
Moderator at 345.15 K with a density of 1.086 g/cc
Coolant at 561.65 K with a density of 0.806 g/cc
Pressure tube at 600 K with density of 6.3 g/cc
Calandria tube at 300 K with density of 6.4 g/cc

Also, your coolant and moderator purity needs to be adjusted. This is very important in CANDU lattice physics. Both are mixtures of light and heavy water. The moderator is 99.925 wt % D2O, while the coolant is 98.39 wt. %D2O. The nominal Boron concentration is 0.1 ppm.

Using these nominal values with the ENDF/B-VII.1 library, you should expect fresh fuel k-eff to be about 1.147 for OpenMC.

@sourena wrote the CNL CANDU model for OpenMC and may be able to provide more insight into the origin of the discrepancy.

1 Like
3

The main issue is with your boundary condition. The CANDU lattice cell has a square boundary of width = 28.575 cm with reflective boundary conditions. You have a cylindrical boundary with vacuum BC. That will probably solve most of the discrepancy.

Depending on how the Serpent model is defined, you may need to refine the materials and dimensions as well, as @AlexandreTrottier pointed out. Double check to make sure they match. The coolant/moderator are the most important, but the pressure tube and calandria tube shouldn’t both be just pure elemental Zr, for example.

Also, for depletion we usually don’t use a linear timestep scale. You might want to start with small, incrementally larger timesteps early on before going to larger steps (something like [0.1, 0.2, 0.5, 1.0, 2.0,…] etc.) Since this is a 2D model, you may want to consider using power density instead of straight up power. For a standard 37-element CANDU bundle, we generally use 31.5 W/g.

Hope that helps.

1 Like
4

Thank you, Dear Expert, for your nice and meaningful answer.

5

Thank you, Dear Expert, for your nice and meaningful answers.

6

Dear expert, I tried to implement what you said in OpenMC. How can I add the Coolant water part to the code section? And when I ran this file, I encountered the following error?

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
# Outer radius of fuel and clad
r_fuel = 0.6122
r_clad = 0.6540
fuel = openmc.Material(name='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

# Update fuel temperature and density
fuel.temperature = 960.0  # in Kelvin
fuel.set_density('g/cm3', 10.52)

# Update heavy water (coolant) temperature and density
heavy_water_coolant = openmc.Material(name='heavy water coolant')
heavy_water_coolant.add_nuclide('H2', 0.799449 )
heavy_water_coolant.add_nuclide('O16', 0.199768 )
heavy_water_coolant.add_nuclide('H1', 0.000783774 )
heavy_water_coolant.set_density('g/cm3', 0.806)
heavy_water_coolant.temperature = 561.65
# Update clad material (calandria tube) temperature and density
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.temperature = 300.0  # in Kelvin
clad_calandria.set_density('g/cm3', 6.4)
#*********************************************
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.temperature = 600.0  # in Kelvin 
clad.set_density('g/cm3', 6.3)
heavy_water = openmc.Material(name='heavy water')
# Update heavy water (moderator) purity
heavy_water.add_nuclide('H2', 0.798895 )
heavy_water.add_nuclide('O16', 0.201016 )
heavy_water.add_nuclide('H1', 8.96e-5 )
# Update heavy water (moderator) temperature and density
heavy_water.temperature = 345.15  # in Kelvin 
heavy_water.set_density('g/cm3', 1.086)
heavy_water.add_s_alpha_beta('c_D_in_D2O')
heavy_water.set_density('g/cm3', 1.1)
# PressMtube and calendria radii
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 )

pressure_tube_ir = 5.16890
pressure_tube_or = 5.60320
calendria_ir = 6.44780
calendria_or = 6.58750

# 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)

fuel_cell = openmc.Cell(fill=fuel, region=-surf_fuel)
clad_cell = openmc.Cell(fill=clad, region=+surf_fuel)

pin_universe = openmc.Universe(cells=(fuel_cell, clad_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
        
        # Create each fuel pin -- note that we explicitly assign an ID so 
        # 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])
geometry = openmc.Geometry(root_universe)
geometry.export_to_xml()

materials = openmc.Materials([fuel, clad, heavy_water])
materials.export_to_xml()
################################################################################
# SIMULATION SETTINGS

settings = openmc.Settings()
settings.particles = 4000000
settings.inactive = 20
settings.batches = 220
settings.statepoint = {'batches': []}
settings.output = {'tallies': False}
settings.temperature = {
    'default': 293.15,
    'method': 'interpolation',
}
settings.export_to_xml()

################################################################################
# create a Model object
model = openmc.Model(geometry=geometry, materials=materials, settings=settings)
# DEPLETION

cumulative_days = np.array([0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0])
timesteps = np.diff(cumulative_days)
# Adjusting power density to 31.5 W/g
power = 31.5  # in W/g

data_path = Path('/mnt/d/openmc/examples/depletion-comparison/data/depletion')
file_path = data_path / 'serpent_fissq.json'

with file_path.open() as fh:
    serpent_q = json.load(fh)


# Create depletion operator
chain_file = data_path / 'chain_simple.xml'
op = openmc.deplete.CoupledOperator(model, chain_file,
    fission_q=serpent_q,
    fission_yield_mode="average"
)

# Execute depletion using integrator
integrator = openmc.deplete.PredictorIntegrator(
    op, timesteps, power, timestep_units='d')
integrator.integrate()

 Reading settings XML file...
 Reading cross sections XML file...
 Reading materials XML file...
 Reading geometry XML file...
 ERROR: Could not find material 6 specified on cell 45

Serpent adding coolant water sample

% --- Fuel pin:



pin 1 

fuel  0.6122  

clad  0.6540  

cool
7

The error you are seeing is because your materials list did not include all of the materials in your model. Each Material you define should be included in that list (clad_calandria, heavy_water_coolant, etc.).

1 Like