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