diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 4a4c7b1f8..f6feb6d24 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -33,12 +33,11 @@ material compositions over time. Each method appears as a different class. For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation using the CE/CM algorithm (deplete over a timestep using the middle-of-step reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to -one of these functions along with the power level and timesteps:: +one of these functions along with the timesteps and power level:: - power = 1200.0e6 - days = 24*60*60 - timesteps = [10.0*days, 10.0*days, 10.0*days] - openmc.deplete.CECMIntegrator(op, power, timesteps).integrate() + power = 1200.0e6 # watts + timesteps = [10.0, 10.0, 10.0] # days + openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate() The coupled transport-depletion problem is executed, and once it is done a ``depletion_results.h5`` file is written. The results can be analyzed using the @@ -67,7 +66,7 @@ the energy deposited during a transport calculation will be lower than expected. This causes the reaction rates to be over-adjusted to hit the user-specific power, or power density, leading to an over-depletion of burnable materials. -There are some remedies. First, the fission Q values can be directly set in a +There are some remedies. First, the fission Q values can be directly set in a variety of ways. This requires knowing what the total fission energy release should be, including indirect components. Some examples are provided below:: @@ -99,11 +98,11 @@ Local Spectra and Repeated Materials ------------------------------------ It is not uncommon to explicitly create a single burnable material across many locations. -From a pure transport perspective, there is nothing wrong with creating a single +From a pure transport perspective, there is nothing wrong with creating a single 3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in every fuel pin in an assembly or even full core problem. This certainly expedites the model making process, but can pose -issues with depletion. -Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using +issues with depletion. +Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using a single set of reaction rates, and produce a single new composition for the next time step. This can be problematic if the same ``fuel_3`` is used in very different regions of the problem. diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index 013a2469e..f0387e066 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -1,25 +1,7 @@ import openmc import openmc.deplete -import numpy as np import matplotlib.pyplot as plt -############################################################################### -# Simulation Input File Parameters -############################################################################### - -# OpenMC simulation parameters -batches = 100 -inactive = 10 -particles = 1000 - -# Depletion simulation 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_simple.xml' -power = 174 # W/cm, for 2D simulations only (use W for 3D) - ############################################################################### # Load previous simulation results ############################################################################### @@ -37,31 +19,34 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5") ############################################################################### # Instantiate a Settings object, set all runtime parameters -settings_file = openmc.Settings() -settings_file.batches = batches -settings_file.inactive = inactive -settings_file.particles = particles +settings = openmc.Settings() +settings.batches = 100 +settings.inactive = 10 +settings.particles = 10000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings_file.source = openmc.source.Source(space=uniform_dist) +settings.source = openmc.source.Source(space=uniform_dist) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] -settings_file.entropy_mesh = entropy_mesh +settings.entropy_mesh = entropy_mesh ############################################################################### # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file, - previous_results) +# Create depletion "operator" +chain_file = './chain_simple.xml' +op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results) # Perform simulation using the predictor algorithm -integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) +time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days +power = 174 # W/cm, for 2D simulations only (use W for 3D) +integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d') integrator.integrate() ############################################################################### @@ -77,27 +62,28 @@ 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)') +# Obtain Xe135 capture reaction rate as a function of time +time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### # Generate plots ############################################################################### +days = 24*60*60 plt.figure() -plt.plot(time/(24*60*60), keff, label="K-effective") +plt.plot(time/days, 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.plot(time/days, 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.plot(time/days, Xe_capture, label="Xe135 capture") plt.xlabel("Time (days)") plt.ylabel("RR (-)") plt.show() diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 9933edd48..ce4b0cf5e 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -1,47 +1,31 @@ +from math import pi + import openmc import openmc.deplete -import numpy as np import matplotlib.pyplot as plt -############################################################################### -# Simulation Input File Parameters -############################################################################### - -# OpenMC simulation parameters -batches = 100 -inactive = 10 -particles = 1000 - -# Depletion simulation 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_simple.xml' -power = 174 # W/cm, for 2D simulations only (use W for 3D) - ############################################################################### # Define materials ############################################################################### # Instantiate some Materials and register the appropriate Nuclides -uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment') +uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') uo2.set_density('g/cm3', 10.29769) uo2.add_element('U', 1., enrichment=2.4) uo2.add_element('O', 2.) -uo2.depletable = True -helium = openmc.Material(material_id=2, name='Helium for gap') +helium = openmc.Material(name='Helium for gap') helium.set_density('g/cm3', 0.001598) helium.add_element('He', 2.4044e-4) -zircaloy = openmc.Material(material_id=3, name='Zircaloy 4') +zircaloy = openmc.Material(name='Zircaloy 4') zircaloy.set_density('g/cm3', 6.55) -zircaloy.add_element('Sn', 0.014 , 'wo') +zircaloy.add_element('Sn', 0.014, 'wo') zircaloy.add_element('Fe', 0.00165, 'wo') -zircaloy.add_element('Cr', 0.001 , 'wo') +zircaloy.add_element('Cr', 0.001, 'wo') zircaloy.add_element('Zr', 0.98335, 'wo') -borated_water = openmc.Material(material_id=4, name='Borated water') +borated_water = openmc.Material(name='Borated water') borated_water.set_density('g/cm3', 0.740582) borated_water.add_element('B', 4.0e-5) borated_water.add_element('H', 5.0e-2) @@ -52,87 +36,62 @@ borated_water.add_s_alpha_beta('c_H_in_H2O') # Create geometry ############################################################################### -# Instantiate ZCylinder surfaces -fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR') -clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR') -clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR') -left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left') -right = openmc.XPlane(surface_id=5, x0=0.62992, name='right') -bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom') -top = openmc.YPlane(surface_id=7, y0=0.62992, name='top') +# Define surfaces +pitch = 1.25984 +fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR') +clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR') +clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR') +box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective') -left.boundary_type = 'reflective' -right.boundary_type = 'reflective' -top.boundary_type = 'reflective' -bottom.boundary_type = 'reflective' +# Define cells +fuel = openmc.Cell(fill=uo2, region=-fuel_or) +gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir) +clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or) +water = openmc.Cell(fill=borated_water, region=+clad_or & box) -# Instantiate Cells -fuel = openmc.Cell(cell_id=1, name='cell 1') -gap = openmc.Cell(cell_id=2, name='cell 2') -clad = openmc.Cell(cell_id=3, name='cell 3') -water = openmc.Cell(cell_id=4, name='cell 4') - -# Use surface half-spaces to define regions -fuel.region = -fuel_or -gap.region = +fuel_or & -clad_ir -clad.region = +clad_ir & -clad_or -water.region = +clad_or & +left & -right & +bottom & -top - -# Register Materials with Cells -fuel.fill = uo2 -gap.fill = helium -clad.fill = zircaloy -water.fill = borated_water - -# Instantiate Universe -root = openmc.Universe(universe_id=0, name='root universe') - -# Register Cells with Universe -root.add_cells([fuel, gap, clad, water]) - -# Instantiate a Geometry, register the root Universe -geometry = openmc.Geometry(root) +# Define overall geometry +geometry = openmc.Geometry([fuel, gap, clad, water]) ############################################################################### # Set volumes of depletable materials ############################################################################### -# Compute cell areas -area = {} -area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2 - -# Set materials volume for depletion. Set to an area for 2D simulations -uo2.volume = area[fuel] +# Set material volume for depletion. For 2D simulations, this should be an area. +uo2.volume = pi * fuel_or.r**2 ############################################################################### # 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 +settings = openmc.Settings() +settings.batches = 100 +settings.inactive = 10 +settings.particles = 1000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings_file.source = openmc.source.Source(space=uniform_dist) +settings.source = openmc.source.Source(space=uniform_dist) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] -settings_file.entropy_mesh = entropy_mesh +settings.entropy_mesh = entropy_mesh ############################################################################### # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file) +# Create depletion "operator" +chain_file = './chain_simple.xml' +op = openmc.deplete.Operator(geometry, settings, chain_file) # Perform simulation using the predictor algorithm -integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) +time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days +power = 174 # W/cm, for 2D simulations only (use W for 3D) +integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d') integrator.integrate() ############################################################################### @@ -148,27 +107,28 @@ 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)') +# Obtain Xe135 capture reaction rate as a function of time +time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### # Generate plots ############################################################################### +days = 24*60*60 plt.figure() -plt.plot(time/(24*60*60), keff, label="K-effective") +plt.plot(time/days, 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.plot(time/days, n_U235, label="U235") plt.xlabel("Time (days)") plt.ylabel("n U5 (-)") plt.show() plt.figure() -plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption") +plt.plot(time/days, Xe_capture, label="Xe135 capture") plt.xlabel("Time (days)") plt.ylabel("RR (-)") plt.show()