Issue with Depletion Calculations

[isak]

Hi! I am having trouble with depletion calculations. I created a typical cubic reactor and I’m trying to perform depletion calculations. Below I have attached some of the important settings and geometry I’ve done; we are running 120 batches. The fuel I’m running the simulations for contains a great percentage Uranium-235 but is contaminated by other myriad of isotopes all over the periodic table. My main issue is that when I first run OpenMC I get a supercritical value (~1.10), but when I perform depletion calculations I get two files with subcritical values (0.62, 0.53: i.e., the numbers are not consistent). And even the graph below shows that our initial keff is past 1.6, not even being consistent with the initial OpenMC run. I have had a graduate student mentor look at it and we cannot pinpoint the issue. Some questions we have asked ourselves are: why is the initial keff changing when we change the volume? Why is the steady-state keff different from the initial burnup core keff, which should be identical? Manipulating our inputs, why is the keff only consistent when the fuel comprises only three nuclides (O16, U235, U238)?
We might think there is an issue with how Jupyter performs the depletion calculations, but we’re not sure. Any help is greatly appreciated and your luck will triple!! We’ve been dealing with this issue for a while with no clear solution. :'^(

The depletion looks like (I would’ve posted a screenshot but this post only allowed one image):

#DEPLETION
import openmc.deplete
model = openmc.Model(geometry=geometry, materials=materials,  settings=settings)
operator = openmc.deplete.CoupledOperator(model, "/Users/kdoll/openmc/lib/chain_endfb71_pwr.xml")

#We will take depletion step sizes of 30days
# and instruct OpenMC to re-run a transport simulation every 30 days until we have modeled the problem over a six month cycle.
#The depletion interface expects the time to be given in seconds, so we will have to convert.
total_power = 3400*10**3 #MWthermal
time_steps = [300] * 1
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, total_power, timestep_units='d')
integrator.integrate()

[wahidluthfi]

Dear Isak, welcome to the openmc community
About the difference in calculated k-eff between the static non-depletion calculation and the initial keff at step-0 of depletion calculation, I think you need to check whether the openmc model generated for depletion calculation differs from the non-depletion calculation. Have you tried to rerun the calculation, either static (non-depletion) or depletion, into separate new folders? because sometimes our working directory is too crowded and confusing.
The material volume could affect our calculations since the total heavy metal mass is correlated with the HM material volume. So for the same amount of power generated, more heavy metal means a longer cycle length since your core can maintain reactivity longer.
About your depletion input, you are using 3400E+3 as total power, 3.4MW for your whole core model right? so your fuel volume should be the total fuel volume in your core model. But if you use a 1/4 core model, then the core power should be 1/4 nominal power with the fuel volume also reduced to 1/4 total fuel volume.
Also, regarding the fuel composition vs calculation consistency, we can’t check it because we didn’t have your model (material composition), but in general, this should be consistent with a small deviation.

[isak]

Hello! Thank you so much for your response. let me show you how I did the geometry.

If this makes a difference in what the issue could be, please let me know!

[isak]

I’m also happy to share the jupyter notebook if necessary!

[wahidluthfi]

Hi Isak, I think it would be great if you wanted to share the notebook. but if you don’t want because it was restricted, I think the fuel composition itself would be good since we can try by calculating kinf with reflective geometry.