PWR UO2 Benchmark

world3 · January 15, 2022, 2:14pm

Hello Paul,
@paulromano I’m sorry to bother you again.Recently,I have used OpenMC to validate a PWR Assembly benchmark problem(Akio YAMAMOTO,Tadashi IKEHARA,Takuya ITO and Etsuro SAJI,Benchmark Problem Suite for Reactor Physics Study of LWR Next Generation Fuels).
There is a little difference with other procedures(I use the ‘chain_casl_pwr.xml’ chain):

Here is my Python Script:

import openmc
import os
from openmc import geometry
import openmc.deplete


dress_work='Verification'
os.chdir(dress_work)

uo2=openmc.Material(1,'UO2 Fuel',900.0)
uo2.add_nuclide('U235',1.5122E-03)
uo2.add_nuclide('U238',2.1477E-02)
uo2.add_nuclide('O16',4.5945E-02)
uo2.set_density('g/cm3',10.3)
uo2.volume=123.7178287866994403
zr=openmc.Material(2,'Structural Material',600.0)
zr.add_nuclide('Zr90',2.2200e-2)
zr.add_nuclide('Zr91',4.8280e-3)
zr.add_nuclide('Zr92',7.3713e-3)
zr.add_nuclide('Zr94',7.5006e-3)
zr.add_nuclide('Zr96',1.2070e-3)
zr.set_density('g/cm3',6.53)

water=openmc.Material(3,'Moderator Material',600.0)
water.add_nuclide('H1',2)
water.add_nuclide('O16',1)
water.add_s_alpha_beta("c_H_in_H2O")
water.set_density ('g/cm3',0.66)

gd=openmc.Material(4,'Gd bearing fuel rod',900.0)
gd.add_nuclide('U235',8.1312E-04)
gd.add_nuclide('U238',1.9268E-02)
gd.add_nuclide('Gd154',7.1289E-05)
gd.add_nuclide('Gd155',4.8938E-04)
gd.add_nuclide('Gd156',6.8028E-04)
gd.add_nuclide('Gd157',5.2077E-04)
gd.add_nuclide('Gd158',8.2650E-04)
gd.add_nuclide('Gd160',7.2761E-04)
gd.add_nuclide('O16',4.5130E-02)
gd.set_density('g/cm3',10)
gd.volume=17.06452810851026763
materilas=openmc.Materials([uo2,zr,water,gd])
materilas.export_to_xml()
# UO2 Pin Cell
uo2_radius=openmc.ZCylinder(r=0.412)
clad_radius_uo2=openmc.ZCylinder(r=0.476)
uo2_cell=openmc.Cell(1,'UO2 Fuel',uo2,-uo2_radius)
clad_cell_uo2=openmc.Cell(2,'UO2 Clad',zr,+uo2_radius&-clad_radius_uo2)
uo2_out_cell=openmc.Cell(3,'UO2 Out Cell',water,+clad_radius_uo2)

u1=openmc.Universe(1,'UO2 Pin',[uo2_cell,clad_cell_uo2,uo2_out_cell])
# Gd bearing Cell
gd_radius=openmc.ZCylinder(r=0.412)
clad_radius_gd=openmc.ZCylinder(r=0.476)
gd_cell=openmc.Cell(4,'Gd Bearing',gd,-gd_radius)
clad_cell_gd=openmc.Cell(5,'Gd clad',zr,+gd_radius&-clad_radius_gd)
gd_out_cell=openmc.Cell(6,'Gd Out Cell',water,+clad_radius_gd)

u2=openmc.Universe(2,'Gd Bearing Pin',[gd_cell,clad_cell_gd,gd_out_cell])
# Thimble
guide_inner=openmc.ZCylinder(r=0.57)
guide_outer=openmc.ZCylinder(r=0.61)
guide_cell=openmc.Cell(7,'Guide Cell',water,-guide_inner)
guide_clad=openmc.Cell(8,'Guide clad',zr,+guide_inner&-guide_outer)
guide_out=openmc.Cell(9,'Guide Out Cell',water,+guide_outer)
u3=openmc.Universe(3,'Guide Universe',[guide_cell,guide_clad,guide_out])
# Assembly
lattice=openmc.RectLattice(4,'Fuel Lattice')
lattice.pitch=(1.265,1.265)
lattice.lower_left=[-1.265*17/2]*2
lattice.universes=[[u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,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,u2,u3,u1,u2,u3,u2,u1,u3,u2,u1,u2,u1,u1],
                   [u1,u1,u1,u3,u1,u1,u1,u1,u1,u1,u1,u1,u1,u3,u1,u1,u1],
                   [u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1],
                   [u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1],
                   [u1,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u1],
                   [u1,u1,u2,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u2,u1,u1],
                   [u1,u1,u3,u1,u1,u3,u2,u1,u3,u1,u2,u3,u1,u1,u3,u1,u1],
                   [u1,u1,u2,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u1,u2,u1,u1],
                   [u1,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u1],
                   [u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1,u3,u1,u1],
                   [u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1,u1,u2,u1,u1],
                   [u1,u1,u1,u3,u1,u1,u1,u1,u1,u1,u1,u1,u1,u3,u1,u1,u1],
                   [u1,u1,u2,u1,u2,u3,u1,u2,u3,u2,u1,u3,u2,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]]
box=openmc.rectangular_prism(21.505,21.505,'z',(0,0),'reflective')
assembly=openmc.Cell(10,'Assembly',lattice,box)
root_universe=openmc.Universe(5,'Root Universe',[assembly])
geometry=openmc.Geometry(root_universe)
geometry.export_to_xml()
point = openmc.stats.Point((0, 0, 0))
source = openmc.Source(space=point)
settings = openmc.Settings()
settings.source = source
settings.particles = 1000
settings.inactive = 50
settings.batches =150
settings.export_to_xml()
openmc.run()
operator = openmc.deplete.Operator(geometry,settings,
     diff_burnable_mats=True)
power=37.5*operator.heavy_metal
depletion=[0.1,4.9,5,5,5,10,20,20]
# time_steps=[]
# for i in range(len(depletion)):
#     time_steps.append(depletion[i]*operator.heavy_metal/power*1e3)
# print(time_steps)
# print(operator.heavy_metal)
integrator = openmc.deplete.PredictorIntegrator(operator, depletion, power_density=37.5,timestep_units='MWd/kg')
integrator.integrate()

In addition,I find that the k-eff has a difference between using the ‘openmc.run()’ function and the step of n0 of depletion.
The value of k-eff using ‘openmc.run()’ function is ‘1.13626 +/- 0.00279’,another one is ‘1.13042 +/- 0.00263’.
Thank you!

paulromano · January 15, 2022, 5:47pm

Please read the user’s manual section about energy deposition in depletion. You need to either:

Make sure that the Q values used for fission in OpenMC match whatever is used by the other codes by passing the fission_q argument to Operator, or
Use normalization_mode='energy-deposition' in Operator.

If you want to go with approach 1, a reasonable set of Q values to use can be found here, but you should really read the manuals for these other codes to find out what they assume for the Q value of fission.

Regarding the difference in keff between calling openmc.run and running depletion, those two values you have are not statistically different. If you run more neutrons, you should find that the values are closer.

world3 · January 16, 2022, 1:35pm

@paulromano Hello Paul,I have used energy_mode=“energy-deposition” with more neutrons and get the result:

Maybe there’s something else I haven’t considered?

paulromano · January 16, 2022, 3:43pm

You’re not running very many neutrons per time step, so you may want to increase that as well in case the uncertainties are too high.

The PredictorIntegrator is a very rudimentary time-integrator and will not give accurate results for large steps, so you may want to consider one of the other integrators, such as CECMIntegrator or CELIIntegrator (note that both of those require two transport solutions per time step, so they will be more computationally expensive).

world3 · January 16, 2022, 4:20pm

@paulromano Thank you! I will try to use the EPCRK4Integrator .

world3 · January 18, 2022, 4:31am

@paulromano Hello Paul.I get the result by using EPCRK4Integrator .

When I compared their averages with the benchmark problem I found some discrepancies.

Is there still something wrong with me?

world3 · January 19, 2022, 8:31am

@paulromano Sorry,I showed wrong result.I think this is the final result.

Pranto · January 19, 2022, 10:53am

@world3, Have you considered splitting the gadolinia fuel pin into equal volumetric zones (at least 10 depletion intra-zones are required to produce converged solution) because burnable absorbers experience strong flux gradients until they are all burned away?

world3 · January 19, 2022, 12:42pm

@Pranto I haven’t considered yet.Thanks for your advice.

paulromano · January 19, 2022, 4:09pm

@world3 How many neutrons are you running for each timestep?

world3 · January 19, 2022, 4:29pm

In the past calculation,I used 20000 neutron for each timestep.

world3 · January 19, 2022, 4:35pm

@paulromano Maybe the problem is as @Pranto said I don’t consider the Gd fuel pin separation.

paulromano · January 19, 2022, 4:41pm

That may be the case, but 20,000 neutrons is also pretty low so you’ll probably want to increase the number of particles per batch to reduce the uncertainty on tallies at each timestep.

world3 · January 19, 2022, 4:58pm

@paulromano May I ask how many neutron in general should I simulate at each timestep？

world3 · January 22, 2022, 7:39am

Hello paul,
There are 223 nuclide in the depletion result file.But I find some nuclides for exmple ‘Ag109_m1’,‘Xe135_m1’,‘Nd145’ are not in the cross section library.
What makes me curious is that in every depletion step,these nuclides can get the number.I’m wondering how OpenMC handles this problem.

paulromano · February 1, 2022, 3:43pm

Regarding your first question, unfortunately there’s no easy answer. It all boils down to how much uncertainty you’re willing to live with. The more neutrons you simulate, the lower the uncertainty on your results will be.

For the second question, OpenMC tracks nuclides that don’t have neutron cross sections via the “depletion chain” file, which defines their half-life, decay modes, etc. These nuclides are part of the burnup matrix that is used to solve the Bateman equations. Perhaps you’re pointing out that some of these nuclides in all likelihood have non-negligible neutron cross sections (Xe135_m1, e.g.), but we can’t make cross sections magically appear out of nowhere – that is a problem that nuclear data evaluators have to solve and include in their future libraries.

world3 · February 1, 2022, 4:13pm

@paulromano Thank you very much for your professional answer！

Topic		Replies	Views
Validation database for OpenMC Validation and Verification	1	522	October 5, 2022
Help regarding assembly validation User Support	2	141	April 6, 2020
On the Deviation of Results in Using OpenMC to Produce Multi group Cross Sections for Openmoc in the Presence of Vacuum Boundaries Validation and Verification	0	245	October 10, 2023
OpenMC-MCNP heating validation Validation and Verification	0	102	May 8, 2024
Source code for openmc.examples User Support	4	618	July 30, 2021

PWR UO2 Benchmark

Related topics