Have a trouble during the depletion of the openmc

[gcx]

Hello：
I am verifying benchmark problem of JAEA by OpenMC. But I have a trouble during the depletion of the openmc. I am verifying the uo2 fuel pin problem and the uo2 fuel assembly problem. I find that the keff is not exactly with the depletion. From 5GWd/t, the keff is like some problem.

this is a uo2 fuel assembly problems

This is my uo2 fuel assembly input file

import matplotlib
import openmc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import openmc.deplete
import math

chain_file="chain_endfb71_pwr.xml"   # 输入燃耗链名称

######################################
import os                            #
from os import listdir               #
old_name = chain_file                # 
new_name = "chain"                   #
os.rename(old_name, new_name)        ######################################
for file_name in listdir():          #     删除xml后缀文件以防干扰计算     #
    if file_name.endswith('.xml'):   ######################################
        os.remove(file_name)         #
for file_name in listdir():          #     删除h5后缀文件以防干扰计算     #
    if file_name.endswith('.h5'):   ######################################
        os.remove(file_name)         #
for file_name in listdir():          #     删除png后缀文件以防干扰计算     #
    if file_name.endswith('.png'):   ######################################
        os.remove(file_name)         #
for file_name in listdir():          #     删除csv后缀文件以防干扰计算     #
    if file_name.endswith('.csv'):   ######################################
        os.remove(file_name)         #
os.rename(new_name, old_name)        #
######################################


# Materials definitions

#uo2_fuel
uo2=openmc.Material(name='uo2')
uo2.add_nuclide('U235', 1.5122E-3)
uo2.add_nuclide('U238', 2.1477E-2)
uo2.add_nuclide('O16', 4.5945E-2)
uo2.set_density('atom/b-cm', 6.893420E-02)
uo2.temperature = 900.0

#gd_fuel
gd=openmc.Material(name='gd')
gd.add_nuclide('U235', 8.1312E-4)
gd.add_nuclide('U238', 1.9268E-2)
gd.add_nuclide('Gd154', 7.1289E-5)
gd.add_nuclide('Gd155', 4.8938E-4)
gd.add_nuclide('Gd156', 6.8028E-4)
gd.add_nuclide('Gd157', 5.2077E-4)
gd.add_nuclide('Gd158', 8.2650E-4)
gd.add_nuclide('Gd160', 7.2761E-4)
gd.add_nuclide('O16', 4.5130E-2)
gd.set_density('atom/b-cm', 6.852695E-02)
gd.temperature = 900.0

#fuel clad
fuel_clad = openmc.Material(name='fuel_clad')
fuel_clad.add_nuclide('Zr90',2.2200E-2)
fuel_clad.add_nuclide('Zr91',4.8280E-3)
fuel_clad.add_nuclide('Zr92',7.3713E-3)
fuel_clad.add_nuclide('Zr94',7.5006E-3)
fuel_clad.add_nuclide('Zr96',1.2070E-3)
fuel_clad.set_density('atom/b-cm', 4.310690E-02)
fuel_clad.temperature = 600.0

#coolant
coolant = openmc.Material(name='coolant')
coolant.add_nuclide('H1', 4.4148E-2)
coolant.add_nuclide('O16', 2.2074E-2)
#coolant.add_s_alpha_beta('c_H_in_H2O')
coolant.set_density('atom/b-cm', 6.622200E-02)
coolant.temperature = 600.0

# Instantiate a Materials collection and export to xml
materials = openmc.Materials([fuel_clad,coolant,uo2,gd])
materials.export_to_xml()
# Geometry definitions

# uo2_fuel_pin

fuel_pin_x0 = openmc.XPlane(x0=-0.6325)
fuel_pin_x1 = openmc.XPlane(x0=0.6325)
fuel_pin_y0 = openmc.YPlane(y0=-0.6325)
fuel_pin_y1 = openmc.YPlane(y0=0.6325)

fuel_or1 = openmc.ZCylinder(r=0.238)
fuel_or2 = openmc.ZCylinder(r=0.336)
fuel_or = openmc.ZCylinder(r=0.412)                                                  
clad_ir = openmc.ZCylinder(r=0.412)
clad_or = openmc.ZCylinder(r=0.476)

fuel_region1 = -fuel_or1
fuel_region2 = +fuel_or1 & -fuel_or2
gd_fuel_region = +fuel_or2 & -fuel_or
fuel_region = -fuel_or
clad_region = +clad_ir & -clad_or
coolant_region = +clad_or #& +fuel_pin_x0 & -fuel_pin_x1 & +fuel_pin_y0 & -fuel_pin_y1

uo2_fuel_cell1  = openmc.Cell(fill=uo2, region=fuel_region1)
uo2_fuel_cell2  = openmc.Cell(fill=uo2, region=fuel_region2)
uo2_fuel_cell  = openmc.Cell(fill=uo2, region=gd_fuel_region)
uo2_fuel_cell  = openmc.Cell(fill=uo2, region=fuel_region)
uo2_fuel_clad_cell  = openmc.Cell(fill=fuel_clad, region=clad_region)
uo2_fuel_coolant_cell = openmc.Cell(fill=coolant, region=coolant_region)

###################################################################################################
radii1 = [0.412, 0.476]          
uo2_surfs = [openmc.ZCylinder(r=r) for r in radii1] 
n_fuel_cell1=3 # uo2划分多少区域
uo2_fuel_pin = openmc.model.pin(uo2_surfs, [uo2, fuel_clad, coolant], subdivisions= {0:n_fuel_cell1}) # ring0 or ring1: 分为多少等区域
uo2_fuel_cells = list(uo2_fuel_pin.cells.values())[:n_fuel_cell1]                #列出从最里面数到第十个区域
for cell in uo2_fuel_cells:
    cell.fill = uo2
############################################################################################################

#gd_fuel_pin
#########################################################################################################################################
#gd棒体积划分
#########################################################################################################################################
#radii2 = [0.412, 0.476]  
gd_surfs = [openmc.ZCylinder(r=r) for r in [0.412, 0.476]] 
n_fuel_cell=3                                                                                   # gd划分多少区域
gd_fuel_pin = openmc.model.pin(gd_surfs, [gd, fuel_clad, coolant], subdivisions= {0:n_fuel_cell}) # ring0 or ring1: 分为多少等区域
gd_fuel_cells = list(gd_fuel_pin.cells.values())[:n_fuel_cell]                                    #列出从最里面数到第十个区域
for cell in gd_fuel_cells:
    cell.fill = gd
#########################################################################################################################################



# G/T :RCC guide thimble
G_T_clad_ir = openmc.ZCylinder(r=0.570)
G_T_clad_or = openmc.ZCylinder(r=0.610)

G_T_region = -G_T_clad_ir
G_T_clad_region = +G_T_clad_ir & -G_T_clad_or
G_T_coolant_region = +G_T_clad_or #& +fuel_pin_x0 & -fuel_pin_x1 & +fuel_pin_y0 & -fuel_pin_y1

G_T_cell  = openmc.Cell(fill=coolant, region=G_T_region)
G_T_clad_cell  = openmc.Cell(fill=fuel_clad, region=G_T_clad_region)
G_T_coolant_cell = openmc.Cell(fill=coolant, region=G_T_coolant_region)

G_T = openmc.Universe(cells=[G_T_cell, G_T_clad_cell, G_T_coolant_cell])

# I/T :Instrumentation thimble
I_T_cell  = openmc.Cell(fill=coolant, region=G_T_region)
I_T_clad_cell  = openmc.Cell(fill=fuel_clad, region=G_T_clad_region)
I_T_coolant_cell = openmc.Cell(fill=coolant, region=G_T_coolant_region)

I_T = openmc.Universe(cells=[I_T_cell, I_T_clad_cell, I_T_coolant_cell])

#uo2_fuel_assembly
uo2_fuel_assembly_lat = openmc.RectLattice(name='uo2_fuel_assembly_lat')
uo2_fuel_assembly_lat.dimension= [17, 17]
uo2_fuel_assembly_lat.lower_left = (-10.7525, -10.7525)
uo2_fuel_assembly_lat.pitch = [1.265,1.265]
C1 = uo2_fuel_pin
C2 = gd_fuel_pin
C3 = G_T
C4 = I_T
uo2_fuel_assembly_lat.universes = [[C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1],
                                   [C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1],
                                   [C1,C1,C2,C1,C2,C3,C1,C2,C3,C2,C1,C3,C2,C1,C2,C1,C1],
                                   [C1,C1,C1,C3,C1,C1,C1,C1,C1,C1,C1,C1,C1,C3,C1,C1,C1],
                                   [C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1],
                                   [C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1],
                                   [C1,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C1],
                                   [C1,C1,C2,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C2,C1,C1],
                                   [C1,C1,C3,C1,C1,C3,C2,C1,C4,C1,C2,C3,C1,C1,C3,C1,C1],#9
                                   [C1,C1,C2,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C2,C1,C1],
                                   [C1,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C1],
                                   [C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1,C3,C1,C1],
                                   [C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1,C1,C2,C1,C1],
                                   [C1,C1,C1,C3,C1,C1,C1,C1,C1,C1,C1,C1,C1,C3,C1,C1,C1],
                                   [C1,C1,C2,C1,C2,C3,C1,C2,C3,C2,C1,C3,C2,C1,C2,C1,C1],
                                   [C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1], 
                                   [C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1,C1]]


min_x = openmc.XPlane(x0=-10.7525, boundary_type='reflective')
max_x = openmc.XPlane(x0=+10.7525, boundary_type='reflective')
min_y = openmc.YPlane(y0=-10.7525, boundary_type='reflective')
max_y = openmc.YPlane(y0=+10.7525, boundary_type='reflective')

uo2_fuel_assembly_cell = openmc.Cell(name='uo2_fuel_assembly_cell')
uo2_fuel_assembly_cell.fill = uo2_fuel_assembly_lat
uo2_fuel_assembly_cell.region = +min_x & -max_x & +min_y & -max_y

uo2_fuel_assembly = openmc.Universe(name='uo2_fuel_assembly')
uo2_fuel_assembly.add_cell(uo2_fuel_assembly_cell)


geom=openmc.Geometry(uo2_fuel_assembly)
geom.export_to_xml()

#tallies
tallies = openmc.Tallies()

# Instantiate a tally Mesh
mesh = openmc.RegularMesh()
mesh.dimension = [17 , 17 ]
mesh.lower_left = [-10.7525, -10.7525]
mesh.upper_right = [+10.7525, +10.7525]

# Instantiate tally Filter
mesh_filter = openmc.MeshFilter(mesh)

# Instantiate the Tally
tally = openmc.Tally(name='mesh tally')
tally.filters = [mesh_filter]
tally.scores = ['fission']

# Add tally to collection
tallies.append(tally)

# OpenMC simulation parameters

batches = 200
inactive = 100
particles = 10000

settings = openmc.Settings()
settings.batches = batches
settings.inactive = inactive
settings.particles = particles
# jia su
settings.output = {'tallies': False}
#n = batches
#settings.statepoint = {'batches': range(100, n + 100, 100)}  #As an example, to write a statepoint file every 100 batches:
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-10.7525, -10.7525, -1e50, 10.7525, 10.7525, 1e50]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.Source(space=uniform_dist)

# Tell OpenMC we want to run in eigenvalue mode
settings.run_mode = 'eigenvalue'

#depl
uo2.volume = 232*math.pi*0.412 ** 2
gd.volume = 32*math.pi*0.412 ** 2


chain = openmc.deplete.Chain.from_xml("./chain_endfb71_pwr.xml")
chain.nuclide_dict

# create a model to tie together geometry, materials, settings, and tallies
model = openmc.Model(geometry=geom, materials = materials, settings=settings, tallies = tallies)
model.export_to_xml()

#statepoint_filename = model.run()
op = openmc.deplete.CoupledOperator(model, "./chain_endfb71_pwr.xml",reduce_chain = False , reduce_chain_level=False)

power = 179*264

#使用python脚本转化一下数据

time_steps_1 = [[2.666666667],[130.6666667/5]*5,[133.3333333/2]*2,[133.3333333/1]*1,[133.3333333/1]*1,[266.6666667/1]*1,[533.3333333/2]*2,[533.3333333/2]*2]
time_steps = [float(x) for item in time_steps_1 for x in item]

integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
integrator.integrate()
settings.export_to_xml()
results = openmc.deplete.Results("depletion_results.h5")


i = 0
while i <= 15 :

 # Load the last statepoint file and keff value
 with openmc.StatePoint('openmc_simulation_n'+str(i)+'.h5') as sp:
    # Get the OpenMC pin power tally data
    mesh_tally = sp.get_tally(name='mesh tally')
    fission_rates = mesh_tally.get_values(scores=['fission'])
 # Reshape array to 2D for plotting
 fission_rates.shape = mesh.dimension

 # Normalize to the average pin power
 fission_rates /= np.mean(fission_rates[fission_rates > 0.])
 fission_rates[fission_rates == 0.] = np.nan

 # Plot the pin powers
 plt.figure()
 plt.imshow(fission_rates, interpolation='none', cmap='jet', origin='lower')
 plt.colorbar()
 plt.title('Pin Powers')
 plt.savefig('Pin Powers'+str(i))
 #plt.show()
 
 np.savetxt('Pin Powers'+str(i)+'.csv', fission_rates, delimiter=',')
 i += 1

[wahidluthfi]

Hi gcx,
Sorry, I am not doing some detailed diagnostic to your input script, but I recommend you check the nuclear data being used for the calculated benchmark data you use as reference, double-check your geometry, depletable material volume, timesteps from the burnup value, power being used for depletion (I think you used average power/1cm height of fuel assembly since your volume was based on 1 cm of fuel height), and do some sensitivity calculations for the time integration method as openmc has various methods under the hood and predictor integrator is the fastest one because of its simplicity. Try CELI, CECM, LEQI, etc, hope it could solve your benchmark issue.
https://docs.openmc.org/en/latest/pythonapi/deplete.html

Also, if the benchmark data for nuclide density vs burnup of some heavy metal isotopes exists, you could compare the nuclide density calculated by Openmc to the benchmark data, because from my experience, the time integration method affects nuclide density on various degrees.

[gcx]

Thank you very much！I had tried all of these, but they leaded to some wrong results. I do not know whether the database is the cause, which is ENDF/B7.3 I used. As for the benchmark of JAEA, it was a benchmark in 2000. Some verification results are based on ENDF/B6.3, JEFF3.2 and so on, which are some old database. As for power, I am verifying a 2D problem. So I used an axial line power(W/cm). As for timesteps, I transformed its units from GWd/t to day by the power density(37.5W/gHM).

Best wish! Hope that you can help me!

[wahidluthfi]

Hi gxc, I usually calculate HM mass before determining the time steps duration. this makes some difference in time steps duration as you have been used before (from the image)

I sent you my Jupyter notebook for your case, by adding heavy metal calculation and some changes in time step definitions

I hope this definition of timesteps could improve your calculations, and you might want to double check your input since I have done some calculations using CELI with detailed burnup steps (sorry, our computer is not that fast, so it takes time), and it seems that your current model had something different than benchmark specs which make the kinf drop and rise more than the benchmark data you sent before.

notebook:
benchmarkdepletionJAEA.ipynb (361.0 KB)

Wahid