Problem with depletion calculation using thermal power value as power in MW.

[Sharif Abu Darda]

I have a problem with the depletion calculation. I have generated the chain file with openmc-make-depletion-chain, below is the code of depletion section I am using. I am simulation a whole reactor core. with over 16000 fuel rods. 

import openmc

import openmc.deplete

import numpy as np

import matplotlib.pyplot as plt

# OpenMC simulation parameters

batches = 60

inactive = 10

particles = 1000

##########################################################################################

#                      Depletion Input Parameters

##########################################################################################

time_step = 1*24*60*60 # s

final_time = 5*24*60*60 # s

time_steps = np.full(final_time // time_step, time_step)

chain_file = './chain_endfb71.xml'

power = 330e6 #reactor thermal power was used

###############################################################################

#                     Set volumes of depletable materials

###############################################################################

# Compute cell areas

area = {}

area[Bundle_A_fuel] = np.pi * fuel_A_fuel.coefficients['r'] ** 2

# Set materials volume for depletion.

st_fuel.volume = area[Bundle_A_fuel] * 200

###############################################################################

#                     Transport calculation settings

###############################################################################

# Instantiate a Settings object, set all runtime parameters, and export to XML

settings_file = openmc.Settings()

settings_file.batches = batches

settings_file.inactive = inactive

settings_file.particles = particles

# Create an initial uniform spatial source distribution over fissionable zones

bounds = [-100, -100, -108.3, 100, 100, 91.7]

uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)

settings_file.source = openmc.source.Source(space=uniform_dist)

###############################################################################

#                   Initialize and run depletion calculation

###############################################################################

op = openmc.deplete.Operator(geometry, settings_file, chain_file)

# Perform simulation using the predictor algorithm

integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)

integrator.integrate()

###############################################################################

#                    Read depletion calculation results

###############################################################################

# Open results file

results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")

# Obtain K_eff as a function of time

time, keff = results.get_eigenvalue()

# Obtain U235 concentration as a function of time

time, n_U235 = results.get_atoms('1', 'U235')

# Obtain Xe135 absorption as a function of time

time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')

###############################################################################

#                            Generate plots

###############################################################################

plt.figure()

plt.plot(time/(24*60*60), keff, label="K-effective")

plt.xlabel("Time (days)")

plt.ylabel("Keff")

plt.show()

plt.figure()

plt.plot(time/(24*60*60), n_U235, label="U 235")

plt.xlabel("Time (days)")

plt.ylabel("n U5 (-)")

plt.show()

plt.figure()

plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")

plt.xlabel("Time (days)")

plt.ylabel("RR (-)")

plt.show()

plt.close('all')

The problem is I see the simulation run for the first time with the proper Keff. But the second time the Keff is very low (1/4) of the first run, and it goes down even more with the next run. Until the 3-4 time, I see an error  ERROR: No fission sites banked on MPI rank 0 

Now, If I set the power value to a lower value say (330)W instead of 330MW. I see the simulation run and the Keff result is acceptable. What am I missing here? 

In addition also before each run, I see warning of a material having negative density like

WARNING: nuclide  Th227  in material  1  is negative (density =  -9.760997752436385e-21  at/barn-cm)

WARNING: nuclide  Th228  in material  1  is negative (density =  -9.802581231893582e-21  at/barn-cm)

WARNING: nuclide  Th232  in material  1  is negative (density =  -1.4292281261795273e-18  at/barn-cm)

WARNING: nuclide  Th233  in material  1  is negative (density =  -1.7339315514675744e-20  at/barn-cm)

WARNING: nuclide  Th234  in material  1  is negative (density =  -2.3299574072632586e-17  at/barn-cm)

WARNING: nuclide  Pa233  in material  1  is negative (density =  -3.853794303347148e-21  at/barn-cm)

WARNING: nuclide  U235  in material  1  is negative (density =  -6.965309863028316e-20  at/barn-cm)

[Andrew Johnson]

Sharif,

It is difficult to determine the exact cause, but I believe a likely candidate is the volumes and/or power normalization. When you are computing the volume of your fuel "st_fuel", you take the area of a bundle and multiply it by 200, assuming that is the number of assemblies in your core. The issue could be that the area computed with "area[Bundle_A_fuel] = np.pi * fuel_A_fuel.coefficients['r'] ** 2" is only the area of a single fuel rod, not of all the fuel rods in a single assembly (usually ~250 rods / PWR assembly). This causes your core to be depleted about much faster than anticipated.

This also leads to the negative densities you are seeing, since the solution method used by OpenMC (CRAM) does not ensure positivity of compositions. The negative values could be consistent with your over depletion as well.

I would recommend ensuring the volume of your burnable materials, "st_fuel" are indeed correct. Regards,

Drew