Depletion compare with restart depletion

Validation and Verification
1

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:

image
image935×452 21.7 KB

(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.

1 Like
2

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:

image
image1190×503 15.8 KB

I guess the variance maybe occurring in this?
image
image799×369 80.1 KB

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

3

@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.

2 Likes
4

@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.

image
image1310×551 33.6 KB

5

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.

1 Like
6

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

image
image806×438 43 KB

7

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

image
image526×519 62.2 KB

8

@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?

9

@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())
11

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.

3 Likes
12

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

image
image786×393 36.8 KB

1 Like