1from math import pi 2 3import openmc 4import openmc.deplete 5import matplotlib.pyplot as plt 6 7############################################################################### 8# Define materials 9############################################################################### 10 11# Instantiate some Materials and register the appropriate Nuclides 12uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') 13uo2.set_density('g/cm3', 10.29769) 14uo2.add_element('U', 1., enrichment=2.4) 15uo2.add_element('O', 2.) 16 17helium = openmc.Material(name='Helium for gap') 18helium.set_density('g/cm3', 0.001598) 19helium.add_element('He', 2.4044e-4) 20 21zircaloy = openmc.Material(name='Zircaloy 4') 22zircaloy.set_density('g/cm3', 6.55) 23zircaloy.add_element('Sn', 0.014, 'wo') 24zircaloy.add_element('Fe', 0.00165, 'wo') 25zircaloy.add_element('Cr', 0.001, 'wo') 26zircaloy.add_element('Zr', 0.98335, 'wo') 27 28borated_water = openmc.Material(name='Borated water') 29borated_water.set_density('g/cm3', 0.740582) 30borated_water.add_element('B', 4.0e-5) 31borated_water.add_element('H', 5.0e-2) 32borated_water.add_element('O', 2.4e-2) 33borated_water.add_s_alpha_beta('c_H_in_H2O') 34 35############################################################################### 36# Create geometry 37############################################################################### 38 39# Define surfaces 40pitch = 1.25984 41fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR') 42clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR') 43clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR') 44box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective') 45 46# Define cells 47fuel = openmc.Cell(fill=uo2, region=-fuel_or) 48gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir) 49clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or) 50water = openmc.Cell(fill=borated_water, region=+clad_or & box) 51 52# Define overall geometry 53geometry = openmc.Geometry([fuel, gap, clad, water]) 54 55############################################################################### 56# Set volumes of depletable materials 57############################################################################### 58 59# Set material volume for depletion. For 2D simulations, this should be an area. 60uo2.volume = pi * fuel_or.r**2 61 62############################################################################### 63# Transport calculation settings 64############################################################################### 65 66# Instantiate a Settings object, set all runtime parameters, and export to XML 67settings = openmc.Settings() 68settings.batches = 100 69settings.inactive = 10 70settings.particles = 1000 71 72# Create an initial uniform spatial source distribution over fissionable zones 73bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] 74uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) 75settings.source = openmc.source.Source(space=uniform_dist) 76 77entropy_mesh = openmc.RegularMesh() 78entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] 79entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] 80entropy_mesh.dimension = [10, 10, 1] 81settings.entropy_mesh = entropy_mesh 82 83############################################################################### 84# Initialize and run depletion calculation 85############################################################################### 86 87# Create depletion "operator" 88chain_file = './chain_simple.xml' 89op = openmc.deplete.Operator(geometry, settings, chain_file) 90 91# Perform simulation using the predictor algorithm 92time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days 93power = 174 # W/cm, for 2D simulations only (use W for 3D) 94integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d') 95integrator.integrate() 96 97############################################################################### 98# Read depletion calculation results 99############################################################################### 100 101# Open results file 102results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5") 103 104# Obtain K_eff as a function of time 105time, keff = results.get_eigenvalue() 106 107# Obtain U235 concentration as a function of time 108time, n_U235 = results.get_atoms('1', 'U235') 109 110# Obtain Xe135 capture reaction rate as a function of time 111time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') 112 113############################################################################### 114# Generate plots 115############################################################################### 116 117days = 24*60*60 118plt.figure() 119plt.plot(time/days, keff, label="K-effective") 120plt.xlabel("Time (days)") 121plt.ylabel("Keff") 122plt.show() 123 124plt.figure() 125plt.plot(time/days, n_U235, label="U235") 126plt.xlabel("Time (days)") 127plt.ylabel("n U5 (-)") 128plt.show() 129 130plt.figure() 131plt.plot(time/days, Xe_capture, label="Xe135 capture") 132plt.xlabel("Time (days)") 133plt.ylabel("RR (-)") 134plt.show() 135plt.close('all') 136