Depletion compare with restart depletion

world3 · December 20, 2021, 7:01am

Hello everyone,
I have a question with the OpenMC restart_depletion.I use the OpenMC exmple ‘run_depletion.py’ and ‘restart_depletion.py’ to simulate depletion,but I find they have a little variance in keff.I choose 100 depletion step,and each step for one day;batches for 100 with 10 inactive and particle for 10000.The keff result are as follows:

(The blue line is the exmple ‘run_depletion.py’ keff result,the green line is the exmple ‘restart_depletion.py’ keff result)
Even if the errors is under 0.08%,I am curious about the difference from.
Thanks in advance,
world3.

world3 · December 20, 2021, 2:26pm

The above curve is my result using the depletion
chain ‘chain_casl_pwr.xml’.I find that the curve trend is difference with the curve using the depletion chain ‘chain_simple.xml’.
Like this:

I guess the variance maybe occurring in this?

And there are a lot of warring like this.But in my model,I use the default temparature at 293K.
Thanks in advance.

paulromano · January 9, 2022, 3:38pm

@world3 There is actually a bug in version 0.12.2 of OpenMC whereby the previous depletion results don’t get set correctly. This has since been fixed so the next release of OpenMC won’t have this problem.

Regarding the warnings about probability tables – that’s an issue with the underlying nuclear data libraries and is pretty common with probability tables that were generated with NJOY. OpenMC will zero out entries in the probability table that have negative values and as far as I know it shouldn’t really affect results.

One last note – chain_simple.xml is really only meant for testing purposes and any results you obtain from using it should not be taken seriously as it only contains a few nuclides.

world3 · January 10, 2022, 1:57am

@paulromano Thanks,Paul.I’m appreciate.By the way,I’m wondering why the curve using ‘chain_casl_pwr.xml’ has a downward trend followed by an upward trend.

paulromano · January 10, 2022, 5:13am

To me that just looks like the result of stochastic uncertainty that is inherent in a Monte Carlo simulation. If you were to run this problem with many more particles, you shouldn’t see as much jitter in the results.

world3 · January 10, 2022, 5:32am

@paulromano This group of keff using ‘chain_casl_pwr.xml’ seems strange.

world3 · January 10, 2022, 5:43am

And I use the 17*17 fuel lattice with different burnable materials.Thanks a lot!

paulromano · January 10, 2022, 1:10pm

@world3 I agree it’s a bit strange that k-eff drops so much after the first timestep and then jumps back up. Are you able to share the Python script you are running so I can investigate?

world3 · January 10, 2022, 1:57pm

@paulromano Here is my scipt.Thanks in advance.

import math
import os
from numpy.core.numeric import Infinity
from numpy.lib.utils import source
from openmc import cell, nuclide, region
from openmc.deplete import results
from openmc.deplete.results_list import ResultsList
import openmc.model
import openmc
import openmc.filter
from openmc import geometry
from openmc.model.funcs import subdivide
from openmc.model.model import Model
import openmc.stats
import openmc.statepoint
import openmc.deplete
import numpy as np
from xml.etree.ElementTree import ElementTree
dress_work = '/home/jason/study'
os.chdir(dress_work)
fuel_lat = 17
ctrl_rods_num_s = 24 
axial_div = 11  
fuel_dict = {"U235": 7.67425E+20, "U238": 2.3686E+22, 'O16': 4.89068E+22,
             'I135': 0, 'Xe135': 0, 'Xe136': 0, 'Cs135': 0, 'Gd157': 0, 'Gd156': 0, 'U234': 0, }
rods_dict = {'Ag107': 2.3523E-02, 'Ag109': 2.1854e-02, 'In113': 3.4291E-04, 'In115': 7.6504E-03,
             'Cd108': 2.4221E-05 + 3.4019E-05,
             'Cd110': 3.3991E-04, 'Cd111': 3.4835E-04, 'Cd112': 6.5669E-04, 'Cd113': 3.3257E-04, 'Cd114': 7.8188E-04,
             'Cd116': 2.0384E-04}
fuel_dict_list = []
for i in range((fuel_lat**2-ctrl_rods_num_s)*axial_div):
    fuel_dict_list.append(fuel_dict)

fuel_density = 7.67425E+20 + 2.3686E+22+4.89068E+22
fuel_density_list = []
for i in range((fuel_lat**2-ctrl_rods_num_s)*axial_div):
    fuel_density_list.append(fuel_density)
fuel_volume = math.pi*0.4096**2*10

def settings():
    point = openmc.stats.Point((0, 0, 0))
    source = openmc.Source(space=point)
    settings = openmc.Settings()
    settings.source = source
    settings.particles = 20000
    settings.inactive = 50
    settings.batches = 150
    return settings
def init_depletion_model(aim_rod_height):
    point = openmc.stats.Point((0, 0, 0))
    source = openmc.Source(space=point)
    settings = openmc.Settings()
    settings.source = source
    settings.particles = 10000
    settings.inactive = 50
    settings.batches = 150

    tallies = openmc.Tallies()
    mesh = openmc.RectilinearMesh(name='3d_lattice tally')
    x_grid = []
    y_grid = []
    z_grid = []
    for i in range(18):
        x_grid.append(-10.71 + 1.26 * i)
    mesh.x_grid = x_grid
    for j in range(18):
        y_grid.append(-10.71 + 1.26 * j)
    mesh.y_grid = y_grid
    for k in range(11):
        z_grid.append(-55 + 10 * k)
    mesh.z_grid = z_grid
    energy_filter = openmc.EnergyFilter([0, 0.625, 20.0e6])

    mesh_filter = openmc.MeshFilter(mesh)
    tally = openmc.Tally(name='mesh tally')
    tally.filters = [mesh_filter, energy_filter]
    tally.scores = ['flux', 'fission', 'nu-fission']
    tallies.append(tally)
    tallies.export_to_xml()

    
    fuel_frg = []
    for i in range((fuel_lat**2-ctrl_rods_num_s)*axial_div):
        fuel_frg.append(openmc.Material(material_id=i+1,
                        name='uo2 fragmentation'))
        for f_nu in fuel_dict_list[i].keys():
            fuel_frg[i].add_nuclide(f_nu, fuel_dict_list[i][f_nu])
        fuel_frg[i].set_density("atom/cm3", fuel_density_list[i])
        fuel_frg[i].volume = fuel_volume
        
    ctrl_rods = openmc.Material(material_id=(
        fuel_lat**2-ctrl_rods_num_s)*axial_div+1, name='control rods')
    for r_nu in rods_dict.keys():
        ctrl_rods.add_nuclide(r_nu, rods_dict[r_nu])
    ctrl_rods.set_density('g/cm3', 10.159)
    water = openmc.Material(material_id=(
        fuel_lat**2-ctrl_rods_num_s)*axial_div+2, name='water')
    water.add_nuclide('H1', 2)
    water.add_nuclide('O16', 1)
    water.set_density('g/cm3', 1.0)
    water.add_s_alpha_beta("c_H_in_H2O")
    materials = openmc.Materials([ctrl_rods, water])
    for i in range((fuel_lat**2-ctrl_rods_num_s)*axial_div):
        materials.append(fuel_frg[i])


    ctrl_rods_height_top = 55
    rods_height = ctrl_rods_height_top - aim_rod_height
    if rods_height % 10 == 0 and aim_rod_height == ctrl_rods_height_top:
        mix_height = -5
        position = 0
    elif rods_height % 10 == 0:
        mix_height = -5
        position = 1
    else:
        mix_height = 5-(rods_height % 10)
        position = 0

    # fuel cell
    fuel_radius = openmc.ZCylinder(r=0.4096)
    pin_cell = openmc.Cell(name='pin_cell', fill=fuel_frg, region=-fuel_radius)
    box_out = openmc.Cell(name='out cell', fill=water, region=+fuel_radius)

    # control rods cell
    ctrl_rods_radius = openmc.ZCylinder(r=0.5605)
    ctrl_cell = openmc.Cell(name='control_rods', fill=ctrl_rods)

    rods_region = -ctrl_rods_radius
    ctrl_cell.region = rods_region
    out_region = ~rods_region
    ctrl_cell_outer = openmc.Cell(
        name='outer cell', fill=water, region=out_region)

    # fuel & rods cell
    h = openmc.ZPlane(z0=mix_height)
    water_part_region = -h
    fuel_part = openmc.Cell(
        name='water part', fill=water, region=water_part_region)

    rods_part_region = +h & -ctrl_rods_radius
    rods_part = openmc.Cell(
        name='rods part', fill=ctrl_rods, region=rods_part_region)
    out_rods_part_region = +h & + ctrl_rods_radius
    out_rods_part = openmc.Cell(
        name='out rods part', fill=water, region=out_rods_part_region)

    # water cell
    water_box = openmc.rectangular_prism(1.26, 1.26)
    water_box_region = water_box
    water_cell = openmc.Cell(
        name='water box', fill=water, region=water_box_region)

    u1 = openmc.Universe(name='fuel_universe',
                         cells=[pin_cell, box_out])
    u2 = openmc.Universe(name='rod_universe', cells=[
        ctrl_cell, ctrl_cell_outer])
    u3 = openmc.Universe(name='mix_universe', cells=[
        fuel_part,  rods_part, out_rods_part])
    u4 = openmc.Universe(name='water universe', cells=[water_cell])

    u_fuel = [
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u4, u1, u1, u4,
            u1, u1, u4, u1, u1, u1, u1, u1],
        [u1, u1, u1, u4, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u4, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u4, u1, u1, u4, u1, u1, u4,
            u1, u1, u4, u1, u1, u4, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u4, u1, u1, u4, u1, u1, u1,
            u1, u1, u4, u1, u1, u4, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u4, u1, u1, u4, u1, u1, u4,
            u1, u1, u4, u1, u1, u4, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u4, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u4, u1, u1, u1],
        [u1, u1, u1, u1, u1, u4, u1, u1, u4,
            u1, u1, u4, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1]]

    u_rods = [
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u2, u1, u1, u2,
            u1, u1, u2, u1, u1, u1, u1, u1],
        [u1, u1, u1, u2, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u2, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u2, u1, u1, u2, u1, u1, u2,
            u1, u1, u2, u1, u1, u2, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u2, u1, u1, u2, u1, u1, u1,
            u1, u1, u2, u1, u1, u2, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u2, u1, u1, u2, u1, u1, u2,
            u1, u1, u2, u1, u1, u2, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u2, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u2, u1, u1, u1],
        [u1, u1, u1, u1, u1, u2, u1, u1, u2,
            u1, u1, u2, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1]]

    u_fuel_with_rods = [
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u3, u1, u1, u3,
            u1, u1, u3, u1, u1, u1, u1, u1],
        [u1, u1, u1, u3, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u3, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u3, u1, u1, u3, u1, u1, u3,
            u1, u1, u3, u1, u1, u3, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u3, u1, u1, u3, u1, u1, u1,
            u1, u1, u3, u1, u1, u3, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u3, u1, u1, u3, u1, u1, u3,
            u1, u1, u3, u1, u1, u3, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u3, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u3, u1, u1, u1],
        [u1, u1, u1, u1, u1, u3, u1, u1, u3,
            u1, u1, u3, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1,
            u1, u1, u1, u1, u1, u1, u1, u1],
        [u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1, u1]]
    z_separation = 11
    z_separation_unit = 10
    position = int(z_separation - 1 -
                   (rods_height // z_separation_unit))+position

    lattice = openmc.RectLattice(name='mix_fuel')
    lattice.pitch = (1.26, 1.26, z_separation_unit)
    lattice.lower_left = [-1.26 * 17. / 2.0] * \
        2 + [-z_separation_unit * z_separation / 2.0]

    lattice_list = [0] * z_separation

    lattice_list[position] = u_fuel_with_rods

    for i in range(0, position):
        lattice_list[i] = u_fuel

    for j in range(position + 1, z_separation):
        lattice_list[j] = u_rods
    lattice.universes = lattice_list

    box1 = openmc.rectangular_prism(
        width=21.42, height=21.42, boundary_type='reflective')
    z_min = openmc.ZPlane(z0=-55)
    z_min.boundary_type = 'reflective'
    z_max = openmc.ZPlane(z0=55)
    z_max.boundary_type = 'reflective'
    z_region = +z_min & -z_max & box1
    root_cell = openmc.Cell(
        name='root cell', fill=lattice, region=z_region)
    root_universe = openmc.Universe(
        name='root universe', cells=[root_cell])
    geometry_depletion = openmc.Geometry(root_universe)
    return geometry_depletion
def depletion(settings):
    restart_geometry = init_depletion_model(0)
    operator = openmc.deplete.Operator(restart_geometry, settings,
     energy_mode='energy-deposition',dilute_initial=0)
    time_steps = [1*24*60*60]*5
    integrator = openmc.deplete.PredictorIntegrator(
        operator, time_steps, power_density=200)
    integrator.integrate()
depletion(settings())

paulromano · January 11, 2022, 5:52pm

In the burnup equations, the future concentration of a given nuclide depends on the production rate and the destruction rate. For Xe135, production comes from fission (indirectly) and the destruction is due to absorption. The way that OpenMC estimates the destruction is by setting up reaction rate tallies for the given nuclides – importantly, these reaction rates depend on the actual concentration of the nuclide present in the model. This leads to a potential problem: if a nuclide is not present at the beginning of a depletion simulation, its reaction rate will be tallied as zero. We get around this by adding a very small concentration of the nuclide at t=0 so that we can at least estimate reaction rates (this is the purpose of the dilute_initial parameter). However, in your case, you’ve forced the Xe135 concentration to zero at t=0 by 1) using dilute_initial=0 and 2) explicitly setting the concentration to zero:

The end result is that the first transport simulation gets a Xe135 absorption rate of zero, and hence when we solve the burnup equations, there is only production of Xe135 but no destruction. This results in too much Xe135 after the first timestep, which causes that huge dip in k-eff that you’re seeing. After the second timestep, we have a better estimate of the production/destruction of Xe135 and it settles to its equilibrium value (hence why k-eff jumps back up to where it should have been).

So, in summary: don’t set dilute_initial=0 and get rid of Xe135 (and other zero-concentration nuclides) in your fuel_dict variable and you should see better results.

world3 · January 12, 2022, 6:03am

@paulromano Thank you for your careful explanation. Now I can get the result correctly.I really appereciate it!

Topic		Replies	Views
Issue with Depletion Calculations User Support	4	67	August 16, 2024
Keff remains unchanged in depletion calculation User Support	8	314	June 14, 2024
Restarting depletion calculation, drastic change in Keff User Support	0	400	October 5, 2022
Keff results significantly different from MCNP Validation and Verification	3	687	January 6, 2020
Have a trouble during the depletion of the openmc User Support	3	144	June 24, 2024

Depletion compare with restart depletion

Related topics