-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathTACOCAT.py
211 lines (177 loc) · 9.93 KB
/
TACOCAT.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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from LoadedTACO.src.Clad_Props import Clad_Props
import LoadedTACO.src.HexDhCal as Geom
import LoadedTACO.src.HegNu as Nu
import LoadedTACO.src.TempBulkCal as TempBulk
import TACOCAT_Read_In_File as TCinput
import LoadedTACO.src.Geometry_Value as Geometry
from scipy.integrate import trapz
from scipy.integrate import quad
from LoadedTACO.src.Fuel_Props import Fuel_props
from LoadedTACO.src.Geometry_Value import Core_Geometry
from LoadedTACO.src.Coolant_Value import Coolant
from LoadedTACO.src.Flux_Profiles import Fluxes
from LoadedTACO.src.NusseltNumber import Nusselt
#Assumptions
#1. The core thermal production is assumed to set after heat deposition
#(i.e. gamma isn't relevant)
#----------------------------------------------------------------------------------#
## Print Logicals for saving plots and data files (0 - save, 1 or higher - do not save)
print_logic = TCinput.print_logic
data_logic = TCinput.data_logic
#----------------------------------------------------------------------------------#
## General Core Information
#Read in parameters
Fuel_Type = TCinput.Fuel_Type
Coolant_Type = TCinput.Coolant
Geometry_Type = TCinput.Geometry
Flux_Type = TCinput.Flux
Hc = TCinput.Hc
HotF = TCinput.HotF
Qth = TCinput.Qth
Tbulkin = TCinput.Tbulkin #Bulk Temperature of coolant at the Inlet - C
Uinlet = TCinput.Uinlet #average inlet velocity in a subchannel - m/s
Cladding = TCinput.Clad_Props
Nusselt_Type=TCinput.Nusselt
#Coolant Parameters
Cp = Coolant[Coolant_Type]["Cp"]
rho = Coolant[Coolant_Type]["rho"]
Nu=Nusselt[Nusselt_Type]["NuCorrelation"]
#Geometry Parameters
z = Geometry.z
NFuel = Geometry.NFuel
A_xs = Core_Geometry[Geometry_Type]["CoolantFlowArea"]
#----------------------------------------------------------------------------------#
## Material Properties
Pr = Coolant[Coolant_Type]["Cp"]*Coolant[Coolant_Type]["nu"]*Coolant[Coolant_Type]["rho"]/Coolant[Coolant_Type]["k"] #Prandtl Number Calculation
#Thermal Conductivity of Cladding - W/m-K @ 300 C
kclad = Clad_Props[Cladding]["kclad"] #(Tbulkin + 273.15)
#----------------------------------------------------------------------------------#
## Core Geometry Calculations
# Geometric Calculations
CVol = Core_Geometry[Geometry_Type]["FaceArea"]*Hc #Volume of the core - m^3
#----------------------------------------------------------------------------------#
## Core Parameter Calculations
# Power Density and Linear Generation Calculations
Qavgp = Qth/CVol #Average Power Density Calculation - W/m^3
qlin = Qth/(Geometry.NFuel*Hc) # Average Rod Linear Energy Generation Rate - W/m
qlinHotF = qlin*HotF # Hottest Channel Linear Energy Generation Rate - W/m
qppco = qlin/(np.pi*Geometry.FoCD) # Average heat flux at rod/coolant interface - W/m^2
# Coolant Calculations
mdot = Uinlet*Coolant[Coolant_Type]["rho"]*Core_Geometry[Geometry_Type]["CoolantFlowArea"] # Mass flow rate for the fluid - kg/s
Re = (Uinlet*Core_Geometry[Geometry_Type]["InnerHydraulicDiameter"]/Coolant[Coolant_Type]["nu"])
Pe = Re*Pr # Peclet Number for Fluid
Nu = Nusselt[Nu_correlation]["NuCorrelation"]
h = Nu*Coolant[Coolant_Type]["k"]/Core_Geometry[Geometry_Type]["InnerHydraulicDiameter"] #Heat Transfer Coefficient for Rod Bundles - W/m^2 - C
#Core Temperature Calculations
#Call axial bulk temperature distribution calculation
# Bulk Temperature of Coolant in Average Channel - C
def BulkTemp(Tbulkin,FluxPro,z,NFuel,qlin,Cp,Uinlet,rho,A_xs):
Tbulk = np.zeros(Geometry.steps) #Initialize Bulk Temperature of Coolant - C
for i in range(0,Geometry.steps):
Tbulk[i] = Tbulkin + (np.trapz(FluxPro[0:i+1],z[0:i+1])*NFuel*qlin)/(Cp*Uinlet*rho*A_xs) #Bulk Temperature of Coolant - C
Tbulk[0] = Tbulkin
return Tbulk
# Bulk Temperature of Coolant in Hottest Channel - C
def HotFBulkTemp(Tbulkin,FluxPro,z,NFuel,qlinHotF,Cp,Uinlet,rho,A_xs):
TbulkHotF = np.zeros(Geometry.steps)
for i in range(0,Geometry.steps):
TbulkHotF[i] = Tbulkin + (np.trapz(FluxPro[0:i+1],z[0:i+1])*NFuel*qlinHotF)/(Cp*Uinlet*rho*A_xs) #Bulk Temperature of Coolant - C
TbulkHotF[0] = Tbulkin
return TbulkHotF
FluxPro = Fluxes[Flux_Type]["Profile"]
Tbulk = np.zeros(Geometry.steps)
Tbulk = BulkTemp(Tbulkin,FluxPro,z,NFuel,qlin,Cp,Uinlet,rho,A_xs)
TbulkHotF = np.zeros(Geometry.steps)
TbulkHotF = HotFBulkTemp(Tbulkin,FluxPro,z,NFuel,qlinHotF,Cp,Uinlet,rho,A_xs)
Tcl = np.zeros(Geometry.steps)
TclHotF = np.zeros(Geometry.steps)
for i in range(0,Geometry.steps):
# Centerline Temperature of Fuel in Average Channel - C
Tcl[i] = Tbulk[i] + ((FluxPro[i]*qlin)/(2*np.pi))*((1/(h*(Geometry.FoCD/2))) + (1/(2*Fuel_props[Fuel_Type]["kfuel"])) + (1/kclad)*np.log(Geometry.FoCD/Geometry.FoD))
# Centerline Temperature of Fuel in Hottest Channel - C
TclHotF[i] = TbulkHotF[i] + ((FluxPro[i]*qlinHotF)/(2*np.pi))*((1/(h*(Geometry.FoCD/2))) + (1/(2*Fuel_props[Fuel_Type]["kfuel"])) + (1/kclad)*np.log(Geometry.FoCD/Geometry.FoD))
Tavg = (Tbulk[0] + Tbulk[Geometry.steps-1])/2
THotFavg = (TbulkHotF[0] + TbulkHotF[Geometry.steps-1])/2
Prmax = Coolant[Coolant_Type]["Cpmax"]*Coolant[Coolant_Type]["numax"]*Coolant[Coolant_Type]["rhomax"]/Coolant[Coolant_Type]["kmax"]#Prandtl Number Calculation
#Max and Averaged quantities with calculated outlet temperature
Pravg = (Pr + Prmax)/2 # Average Prandtl Number in a inner channel
rhoavg = (Coolant[Coolant_Type]["rho"] + Coolant[Coolant_Type]["rhomax"])/2 # Average density number in a inner channel
Uoutlet = mdot/(Coolant[Coolant_Type]["rhomax"]*Core_Geometry[Geometry_Type]["CoolantFlowArea"]) # Outlet Velocity in a inner channel
Uavg = (Uinlet + Uoutlet)/2 # Average velocity in a inner channel
Pemax = Prmax*Uoutlet*Geometry.FoCD/Coolant[Coolant_Type]["numax"] # Max Peclet number in a inner channel
Peavg = (Pe + Pemax)/2 # Average Peclect number in a inner channel
Re1 = Peavg/Pravg # Average Reynolds number in a inner channel
if Geometry_Type == "Wire_Wraped_Hexagonal":
M = ((1.034/(Geometry.PtoD**0.124))+(29.7*(Geometry.PtoD**6.94)*Re1**0.086)/(Geom.LeadW(FoCD,WoD)/Geometry.FoCD)**2.239)**0.885 #modifier
else:
M = 1
fsm = 0.316*Re1**(-0.25)
dP = M*fsm*(Hc/Core_Geometry[Geometry_Type]["InnerHydraulicDiameter"])*0.5*rhoavg*Uavg**2
#Heat Load Calculation
QPri = Uavg*rhoavg*Core_Geometry[Geometry_Type]["CoolantFlowArea"]*((Coolant[Coolant_Type]["Cp"] + Coolant[Coolant_Type]["Cpmax"])/2)*(Tbulk[Geometry.steps-1]-Tbulk[0])
#----------------------------------------------------------------------------------#
## Report out - Core Parameters
## Print out
print('--------------------------------------------------------')
print('Heat Load:',QPri/10**6,'MW')
print('Average Power Density:',Qavgp/10**6,'MW/m^3')
print('Average Linear Energy Gen Rate per Rod:',qlin/10**3,'kW/m')
print('Average Heat Flux at Interface:',qppco/10**6,'MW/m^2')
print('Hot Channel Factor:',HotF)
print('Mass Flow Rate:',mdot,'kg/s')
print('Velocity:',Uinlet,'m/s')
print('Inlet Bulk Temperature:',Tbulk[0],'C')
print('Outlet Bulk Temperature:',Tbulk[Geometry.steps-1],'C')
print('Average Bulk Temperature:',Tavg,'C')
print('Outlet Bulk Temperature - Hot Channel',TbulkHotF[Geometry.steps-1],'C')
print('Average Coolant Temperature - Hot Channel:',THotFavg,'C')
print('Highest Centerline Temperature:',Tcl[Geometry.steps-1],'C')
print('Highest Centerline Temperature - Hottest Channel:',TclHotF[Geometry.steps-1],'C')
print('Pressure Drop - Bundle:',dP,'Pa')
print(' Cladding Thermal Conductivity:', kclad, 'W/m-K')
print('--------------------------------------------------------')
#----------------------------------------------------------------------------------#
## Report out - Core Parameters
## Create data files for each run and print out the data
df1 = pd.DataFrame([[Tbulk[0]], [Tbulk[Geometry.steps-1]], [Tavg], [TbulkHotF[Geometry.steps-1]], [THotFavg], [Tcl[Geometry.steps-1]], [TclHotF[Geometry.steps-1]]], index=['Inlet Bulk Temperature', 'Outlet Bulk Temperature', 'Average Bulk Temperature', 'Outlet Bulk Temperature - Hot Channel', 'Average Coolant Temperature - Hot Channel', 'Highest Fuel Centerline Temperature', 'Highest Fuel Centerline Temperature - Hottest Channel'], columns=['Temperature - C'])
df2 = pd.DataFrame({'z': Geometry.z, 'Tbulk': Tbulk, 'TbulkHotF': TbulkHotF, 'Tcl': Tcl, 'TclHotF': TclHotF, 'Flux Profile': FluxPro, 'Fuel Melting Temperature': Fuel_props[Fuel_Type]["Tmelt"], 'Coolant Boiling Temperature': Coolant[Coolant_Type]["Tboil"]})
title_data = TCinput.Reactor_Title + "_table"
if data_logic == 0:
df1.to_excel(title_data + ".xlsx")
df2.to_csv(title_data + '.csv')
#----------------------------------------------------------------------------------#
## Report out - Core Parameters
## Plots
h = 10
w = 8
lw = 2.5
fs = 14
k = 1
plt.figure(k, figsize=(h,w))
plt.plot(Geometry.z,Tbulk,'r--',linewidth = lw,label=r'$T_{bulk}$')
plt.plot(Geometry.z,TbulkHotF,'k--',linewidth = lw, label=r'$T_{bulk-HF}$')
plt.plot(Geometry.z,Tcl, 'g-',linewidth = lw, label=r'$T_{cl}$')
plt.plot(Geometry.z,TclHotF, 'c-',linewidth = lw, label=r'$T_{cl-HF}$')
plt.axhline(y=Coolant[Coolant_Type]["Tboil"], xmin = 0, xmax = 1, color = 'r',linewidth = lw, label='Coolant Boiling Temp')
plt.axhline(y=Coolant[Coolant_Type]["TmeltCoolant"], xmin = 0, xmax = 1, color = 'b',linewidth = lw, label='Coolant Melting Temp')
plt.legend(loc='center left')
plt.xlabel('Axial Height - m',fontsize = fs)
plt.ylabel('Temperature - C',fontsize = fs)
plt.grid()
fig = TCinput.Reactor_Title + '_Axial_Temperatures'
if print_logic == 0:
plt.savefig(fig + '.png', dpi = 300, format = "png",bbox_inches="tight")
plt.savefig(fig + '.eps', dpi = 300, format = "eps",bbox_inches="tight")
plt.savefig(fig + '.svg', dpi = 300, format = "svg",bbox_inches="tight")
k = k + 1
plt.figure(k, figsize=(h,w))
plt.plot(Geometry.z,FluxPro, 'k-')
plt.suptitle('Axial Core Flux Profile')
plt.ylabel('Normalized Height')
plt.xlabel('Normalized Flux')
plt.grid()
k = k + 1
plt.show()