forked from cnerg/OpenMCActivationStudy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SS_Activation_OpenMC.py
170 lines (140 loc) · 5.67 KB
/
SS_Activation_OpenMC.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 8 08:16:12 2024
@author: Anupama Rajendra
"""
import openmc
import openmc.deplete
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import pandas as pd
import random
# Importing Vitamin-J energy group structure:
# This excel file contains the energy bounds of the Vitamin J structure
# Vit_J = pd.read_excel('VitaminJEnergyGroupStructure.xlsx')
# ebounds = Vit_J.iloc[:, 1]
openmc.config['chain_file'] = 'chain_endfb71_sfr.xml'
# Create materials & export to XML:
#Simulating tungsten shell:
W = openmc.Material(material_id=1, name='W_Shell')
W.set_density('g/cm3', 19.28)
W.add_element('W', 1.00)
materials = openmc.Materials([W])
materials.export_to_xml()
# Create geometry
#Spherical shell:
R_1 = 1000
S_1= openmc.Sphere(r=R_1) #sphere of radius 1000cm
inside_sphere_1 = -S_1
outside_sphere_1 = +S_1
R_2 = 1005
S_2 = openmc.Sphere(r=R_2, boundary_type='vacuum')
inside_sphere_2 = -S_2
outside_sphere_2 = +S_2
S_3 = outside_sphere_1 & inside_sphere_2 #filled with tungsten
# Mapping materials to geometry:
Void = openmc.Cell(fill=None, region = inside_sphere_1)
Shell = openmc.Cell(fill=W, region=S_3)
Cells = [Void, Shell]
geometry = openmc.Geometry(Cells)
geometry.export_to_xml()
# Source distribution:
PointSource = openmc.stats.Point(xyz=(0.0, 0.0, 0.0))
Prob = openmc.stats.Discrete(14E+06, 1.0)
# Assign simulation settings
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 1
settings.particles = 10000
settings.source = openmc.Source(space=PointSource, energy=Prob, strength = 1.0, particle = 'neutron')
settings.run_mode = 'fixed source'
settings.export_to_xml()
# Define tallies
neutron_tally = openmc.Tally(name="Neutron tally")
neutron_tally.scores = ['flux', 'elastic', 'absorption']
# Implementing filter for neutron tally through W shell
cell_filter = openmc.CellFilter([Shell])
neutron_tally.filters = [cell_filter]
# Creating a tally to get the flux energy spectrum.
# An energy filter is created to assign to the flux tally using the Vitamin J structure.
energy_filter_flux = openmc.EnergyFilter.from_group_structure("VITAMIN-J-175")
spectrum_tally = openmc.Tally(name="Flux spectrum")
# Implementing energy and cell filters for flux spectrum tally
spectrum_tally.filters = [energy_filter_flux, cell_filter]
spectrum_tally.scores = ['flux']
# Collecting and exporting tallies to .xml
tallies = openmc.Tallies([neutron_tally, spectrum_tally])
tallies.export_to_xml()
model = openmc.model.Model(geometry=geometry,settings=settings)
#Depletion calculation
W.depletable = True
W.volume = 4.0/3.0 * np.pi * (R_2**3 - R_1**3) #volume of W wall material
fluxes, micros = openmc.deplete.get_microxs_and_flux(model, Cells)
operator = openmc.deplete.IndependentOperator(materials, fluxes[0:1], micros[0:1],normalization_mode='source-rate')
# operator = openmc.deplete.CoupledOperator(model, normalization_mode='source-rate')
time_steps = [3e8, 86400, 2.6e6]
source_rates = [1E+18, 0, 0]
integrator = openmc.deplete.PredictorIntegrator(operator=operator, timesteps=time_steps, source_rates=source_rates, timestep_units='s')
integrator.integrate()
#Opening statepoint file to read tallies:
with openmc.StatePoint('statepoint.10.h5') as sp:
fl = sp.get_tally(name="Flux spectrum")
nt = sp.get_tally(name="Neutron tally")
# Get the neutron energies from the energy filter
energy_filter_fl = fl.filters[0]
energies_fl = energy_filter_fl.bins[:, 0]
# Get the neutron flux values
flux = fl.get_values(value='mean').ravel()
#Neutron flux/elastic/absorption tallies:
tal = nt.get_values(value='mean').ravel()
print(tal)
Flux_Data = np.c_[energies_fl, flux]
#Creating an excel file that stores flux data for each energy bin (used as input for ALARA)
#Dividing by volume to obtain proper units of flux (#/cm^2-s)
Flux_Data = np.c_[energies_fl, flux/W.volume]
#ALARA flux inputs go from high energy to low energy
Flux_Data_ALARA = Flux_Data[::-1]
FD_CSV = pd.DataFrame(Flux_Data_ALARA, columns=['Energy [eV]', 'Flux [n-cm/sp]'])
FD_CSV.to_csv('Neutron_Flux.csv', index=False)
Tallies_CSV = pd.DataFrame(tal)
#Creating a csv file that stores total tally value data
Tallies_CSV.to_csv('Tally_Values.csv', index=False)
# Depletion results file
results = openmc.deplete.Results(filename='depletion_results.h5')
# Stable W nuclides present at beginning of operation (will not be plotted)
stable_nuc = ['W180', 'W182', 'W183', 'W184', 'W186']
#Store list of nuclides from last timestep as a Materials object
materials_last = results.export_to_materials(-1)
# Storing depletion data from 1st material
mat_dep = materials_last[0]
# Obtaining the list of nuclides from the results file
nuc_last = mat_dep.get_nuclides()
# Removing stable W nuclides from list so that they do not appear in the plot
for j in stable_nuc :
nuc_last.remove(j)
print(nuc_last)
colors = list(mcolors.CSS4_COLORS.keys())
num_dens= {}
pair_list = {}
with open(r'Densities_CSV.csv', 'a') as density_file:
for nuclide in nuc_last:
plot_color = random.choice(colors)
time, num_dens[nuclide] = results.get_atoms('1', nuclide, nuc_units = 'atom/cm3')
print(time, num_dens[nuclide])
density_file.write(f'{nuclide},')
density_file.write(','.join(map(str, num_dens[nuclide])) + '\n')
plt.plot(time, num_dens[nuclide], marker='.', linestyle='solid', color=plot_color, label=nuclide)
# Adding labels and title
plt.xlabel('Time after beginning of operation [s]')
plt.xlim(1, sum(time_steps)
#plt.gca().set_ylim(bottom=0)
plt.ylabel('Nuclide density [atoms/cm^3]')
plt.xscale("log")
plt.yscale("log")
plt.title('Plot of number density vs time')
# Adding a legend
plt.legend()
plt.savefig('Nuclide_density_OpenMC')
# Display the plot
plt.show()