-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_deplete.py
123 lines (94 loc) · 3.34 KB
/
read_deplete.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
# -*- coding: utf-8 -*-
"""
Created on Sat Apr 15 21:50:33 2023
@author: b9801
"""
import h5py
from pwr1 import get_pins
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 300
import matplotlib.pyplot as plt
import numpy as np
import itertools
colorlist=['#648FFF','#785EF0','#DC267F','#FE6100','#FFB000']
colors = itertools.cycle(colorlist)
next(colors)
ks=list()
fis=['Case A (Reference)','Case D (BL324)','Case E (BL324 with ref. clad thickness)','Case K at +20% power (BL324, 5.95 wt%)']
for fi in ["264","324","324C",'324P']:
f=h5py.File("deplete_"+fi+"/depletion_results.h5")
ks.append([f['eigenvalues'][j][0][0] for j in range(45)])
if fi=="264":
t=[f['time'][j][0]/3600/24/365.25 for j in range(45)]
elif fi != '324P':
dk=[1e5*(k-k0)/k/k0 for k,k0 in zip(ks[-1],ks[0])]
plt.plot(t,dk,color=next(colors),label='BL'+fi)
plt.xlabel('time (yrs)')
plt.ylabel('reactivity difference from reference (pcm)')
plt.legend(loc='lower left')
plt.show()
colors = itertools.cycle(colorlist)
for j in range(4):
plt.plot(t,ks[j],color=next(colors),label=fis[j])
plt.plot([0,7],[1.03,1.03],color='black',linestyle='dashed',label='k-inf=1.03')
plt.xlabel('time (yrs)')
plt.ylabel('k-inf')
plt.legend(loc='upper right')
plt.show()
power=3411e6/193/4/366.0 #Power per cm of the assembly
base_V = (264/4)*0.4095**2*np.pi #cm2 of fuel in assembly
C_mod = ((0.475*np.sqrt(264/324)+0.4095-0.475)/(0.4095*np.sqrt(264/324)))**2 #modifier for the retained clad thickness
HM_rho_5 = 10.4*(235*0.0501+238*(1-0.0501))/(235*0.0501+238*(1-0.0501)+16*2) #HM density for 5% enrichment
HM_rho_6 = 10.4*(235*0.0602+238*(1-0.0602))/(235*0.0602+238*(1-0.0602)+16*2) #HM density for 6% enrichment
base_M = base_V*HM_rho_5
P_mod = HM_rho_6/HM_rho_5
mass=np.array([base_M,base_M,base_M*C_mod,base_M*P_mod])
bu=np.zeros((4,45))
power=[power,power,power,1.2*power]
time=np.array(t)*365.25 #time in days
for j in range(4):
bu[j,:]=power[j]*time[:]/mass[j]/1e3 #1e3 = g to t and W to GW
colors = itertools.cycle(colorlist)
for j in range(4):
plt.plot(bu[j],ks[j],color=next(colors),label=fis[j])
plt.plot([0,100],[1.03,1.03],color='black',linestyle='dashed',label='k-inf=1.03')
plt.xlabel('Burnup (GWd/t))')
plt.ylabel('k-inf')
plt.legend(loc='upper right')
plt.show()
#power
colors = itertools.cycle(colorlist)
p0 = 1.26
p2 = 17/19*p0
P_max=dict()
for n,fi in enumerate(["264","324","324C",'324P']):
P_max[fi]=[]
for j in range(45):
f=h5py.File("deplete_"+fi+"/openmc_simulation_n{}.h5".format(j))
npins = 264 if fi=="264" else 324
ff=np.sqrt(264/npins)
pins = get_pins(p0,p2,ff,npins)
P=[]
total=0
P_total=0
for i in range(len(pins)):
val = f['tallies']['tally {}'.format(i)]['results'][0][0][0]
P_total += val
if pins[i][0] == 0 or pins[i][1] == 0:
val *= 2
total += 0.5
else:
total += 1
P.append(val)
P=np.array(P)
P /= P_total
P *= total
P_max[fi].append(np.max(P))
plt.plot(t,P_max[fi],color=next(colors),label=fis[n])
plt.xlabel('time (yrs)')
plt.ylabel('Assembly power peaking')
plt.legend(loc='upper right')
plt.show()
ks=np.array(ks)
for j in range(4):
print(np.interp(-1.03,-ks[j],t))