我正在尝试运行这个循环；但是，我在代码的第二部分中收到无属性错误。下面是完整的代码（抱歉有点长）。当我运行第一个案例（PWR）时，代码按预期正常执行。但是，当我运行第二种情况（BWR）时，我收到错误，即使它与情况一中的语句完全相同。对此有任何修复或解释吗？谢谢你。
import numpy as np
import math
from iapws import IAPWS97
import matplotlib.pyplot as plt
case = int(input('Which case [1 (PWR) or 2 (BWR)]? '))
if case == 1: # PWR
H = 3.8 # m
He = 3.8 # m
Pitch = 1.25 * 10 ** (-2) # m
Gap_t = 0.00006 # m
D_fuel = 0.0082 # m
k_gap = 0.25 # W/m-K
k_c = 21.5 # W/m-K
k_fuel = 3.6 # W/m-K
T0 = float(278 + 273.15) # K
q0_prime = float(330 * 10 ** (2)) # W/m
P0 = 15 # MPa
MF = float(3460) # kg/m^2-s
D_rod = .0095 # m
R_rod = D_rod / 2
R_fuel = D_fuel / 2
R_gap = R_fuel + Gap_t
R_clad = R_rod
Clad_t = D_rod - D_fuel - Gap_t # m
h0_enthalpy = (IAPWS97(T=T0, P=P0).h) * 10 ** (3)
T_sat0 = IAPWS97(P=P0, x=0).T
g = 9.81 # m/s
# geometry properties
heated_p = math.pi * D_rod
wetted_p = math.pi * D_rod
A_f = (Pitch ** 2) - ((1 / 4) * math.pi * (D_rod ** 2))
D_H = (4 * A_f) / heated_p
# grid setup
grid_points = 100
dz = H / grid_points
z_array = np.arange(0, H, dz)
z_arrayplots = np.arange(0, H, dz)
q_HeatFluxList = []
# defining array of q'' values in list
for z in z_array:
heat_fluxA = (q0_prime / (math.pi * D_rod)) * math.sin(math.pi * (z / He))
q_HeatFluxList.append(heat_fluxA)
q_heat_flux = np.array(q_HeatFluxList)
q_prime = np.zeros(len(z_array))
for i in range(0, len(z_array)):
q_prime[i] = q0_prime * math.sin((np.pi * z_array[i]) / He)
# defining array of h values
h_enthalpy_list = []
h_enthalpy_prefactor = ((heated_p * q0_prime * H) / (A_f * MF * (math.pi ** 2) * D_rod))
for z in z_array:
h_enthalpy = (-h_enthalpy_prefactor * math.cos(math.pi * (z / He))) + h_enthalpy_prefactor + h0_enthalpy
h_enthalpy_list.append(h_enthalpy)
h_enthalpy_array_J = np.array(h_enthalpy_list)
h_enthalpy_array = h_enthalpy_array_J * 10 ** (-3)
P_array = np.zeros(len(z_array))
P_array[0] = P0
T_sat = np.zeros(len(z_array))
T_sat[0] = T_sat0
T_f_array = np.zeros(len(z_array))
T_f_array[0] = T0
Re = np.zeros(len(z_array))
Re_f = np.zeros(len(z_array))
Pr = np.zeros(len(z_array))
k_fluid = np.zeros(len(z_array))
x_array = np.zeros(len(z_array))
xe_array = np.zeros(len(z_array))
frictional = np.zeros(len(z_array))
gravitational = np.zeros(len(z_array))
compressibility = np.zeros(len(z_array))
# Pressure Loop PWR
dp = 0.001
for i in range(0, len(z_array) - 1):
rho_f = IAPWS97(P=P_array[i], x=0).rho
vf = IAPWS97(P=P_array[i], x=0).v
vg = IAPWS97(P=P_array[i], x=1).v
hf_enthalpy = IAPWS97(P=P_array[i], x=0).h
hg_enthalpy = IAPWS97(P=P_array[i], x=1).h
muf = (IAPWS97(P=P_array[i], x=0).mu) * 10 ** (-6)
mug = (IAPWS97(P=P_array[i], x=1).mu) * 10 ** (-6)
k_fluid[i] = IAPWS97(P=P_array[i], T=T_f_array[i]).k
Pr[i] = IAPWS97(P=P_array[i], h=h_enthalpy_array[i]).Liquid.Prandt
x_array[i] = 0
xe_array[i] = (h_enthalpy_array[i] - hf_enthalpy) / (hg_enthalpy - hf_enthalpy)
rho_m = 1 / ((x_array[i] * vg) + ((1 - x_array[i]) * vf))
mu_m = 1 / ((x_array[i] / mug) + ((1 - x_array[i]) / muf))
Re[i] = (MF * D_H) / (mu_m * 10 ** 6) # convert mu to Pa/s
f = 0.079 * (Re[i] ** -0.25) * (mu_m / muf)
Tau = (1 / 2) * f * ((MF ** 2) / rho_m)
Re_f[i] = Re[i]
vf_plus_dP = IAPWS97(P=P_array[i] + dp, x=0).v
vf_minus_dP = IAPWS97(P=P_array[i] - dp, x=0).v
ddP_vf = (vf_plus_dP - vf_minus_dP) / (2 * (dp * 10 ** 6))
frictional[i] = (Tau * wetted_p) / A_f
gravitational[i] = g * rho_f
compressibility[i] = (MF ** 2) * (ddP_vf)
dPdz_num = (frictional[i] + gravitational[i]) # Pa/m
dPdz_denom = 1 + compressibility[i] # Pa/m
dPdz = -dPdz_num / dPdz_denom # Pa/m
P_array[i + 1] = P_array[i] + ((dPdz * dz) * 10 ** (-6))
T_f_array[i + 1] = IAPWS97(P=P_array[i + 1], h=h_enthalpy_array[i + 1]).T
T_sat[i + 1] = IAPWS97(P=P_array[i + 1], x=0).T
# final calc for final value of quality and void fraction because loop stops before these
hf_final = IAPWS97(P=P_array[-1], x=0).h
hg_final = IAPWS97(P=P_array[-1], x=1).h
muf_final = (IAPWS97(P=P_array[-1], x=0).mu) * 10 ** (-6)
mug_final = (IAPWS97(P=P_array[-1], x=1).mu) * 10 ** (-6)
k_fluid[-1] = IAPWS97(P=P_array[-1], T=T_f_array[-1]).k
xe_array[-1] = (h_enthalpy_array[-1] - hf_final) / (hg_final - hf_final)
# fuel and clad temps
T_C_Outer = np.zeros(len(z_array))
mu_m_final = 1 / ((x_array[-1] / mug_final) + ((1 - x_array[-1]) / muf_final))
Re_f[-1] = (MF * D_H) / (muf_final * 10 ** 6)
Pr[-1] = IAPWS97(P=P_array[i], h=h_enthalpy_array[i]).Liquid.Prandt
h_HT = 0.023 * (Re_f[0] ** 0.8) * (Pr[0] ** 0.4) * (k_fluid[0] / D_H)
T_C_Outer[0] = (q_heat_flux[0] + (h_HT * T_f_array[0])) / h_HT
for i in range(0, len(z_array) - 1):
h_HT = 0.023 * (Re_f[i + 1] ** 0.8) * (Pr[i + 1] ** 0.4) * (k_fluid[i + 1] / D_H)
T_C_Outer[i + 1] = (q_heat_flux[i + 1] + (h_HT * T_f_array[i + 1])) / h_HT
q_triple_prime = (q_prime * 4) / (np.pi * (D_fuel ** 2))
T_C_Inner = np.zeros(len(z_array))
T_F_Outer = np.zeros(len(z_array))
T_F_Center = np.zeros(len(z_array))
for i in range(0, len(z_array)):
C1 = -((q0_prime * R_clad) / (k_c * heated_p)) * np.sin(np.pi * (z_array[i] / H))
C2 = T_C_Outer[i] - (C1 * np.log(R_clad))
T_C_Inner[i] = (C1 * np.log(R_gap)) + C2
C3 = (k_c / k_gap) * C1
C4 = T_C_Inner[i] - (C3 * np.log(R_gap))
T_F_Outer[i] = (C3 * np.log(R_fuel)) + C4
C6 = T_F_Outer[i] + ((q_triple_prime[i] * (R_fuel ** 2)) / (4 * k_fuel))
T_F_Center[i] = C6
CL_max = np.amax(T_F_Center)
index = np.where(T_F_Center == CL_max)
z_CL_max = z_array[index]
Clad_max = np.amax(T_C_Inner)
index = np.where(T_C_Inner == Clad_max)
z_Clad_max = z_array[index]
plt.figure(1)
plt.plot(T_C_Outer, z_arrayplots, label='Clad Outer Surface Temp')
plt.plot(T_C_Inner, z_arrayplots, label='Clad Inner Surface Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempClad.png", dpi=600)
plt.figure(2)
plt.plot(T_C_Outer, z_arrayplots, label='Clad Outer Surface Temp')
plt.plot(T_C_Inner, z_arrayplots, label='Clad Inner Surface Temp')
plt.plot(T_F_Outer, z_arrayplots, label='Fuel Outer Surface Temp')
plt.plot(T_F_Center, z_arrayplots, label='Fuel Centerline Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempFuelAndClad.png", dpi=600)
# radial calcs
T_array_A = [T_F_Center[25], T_F_Outer[25], T_C_Inner[25], T_C_Outer[25]]
T_array_B = [T_F_Center[49], T_F_Outer[49], T_C_Inner[49], T_C_Outer[49]]
T_array_C = [T_F_Center[53], T_F_Outer[53], T_C_Inner[53], T_C_Outer[53]]
r_array = [0, R_fuel, R_gap, R_clad]
plt.figure(3)
plt.plot(r_array, T_array_A, label='z = -H/4 = -0.9 m')
plt.plot(r_array, T_array_B, label='z = 0 m')
plt.plot(r_array, T_array_C, '--', label='z = zmax = 0.108 m')
plt.legend(loc='upper left')
plt.ylabel("Temperature [K]")
plt.xlabel("Radius r [m]")
plt.savefig("TempRadial.png", dpi=600)
# critical heat flux and DNBR
P_array_DNBR = np.delete(P_array, 0)
q_heat_flux_DNBR = np.delete(q_heat_flux, 0)
z_arrayplots_DNBR = np.delete(z_arrayplots, 0)
G_Mlbs = MF * (((2.20462 * 10 ** (-6)) * 3600) / 10.7639)
q_heat_flux_MBtu = q_heat_flux[1:] * 3.41 * (1 / 1000000) * (1 / 10.7639)
P_c = 22.064 # https://nuclearstreet.com/nuclear-power-plants/w/nuclear_power_plants/features-of-pressurized-water-reactors
P_crit = P_array_DNBR / P_c
P1 = 0.5328
P2 = 0.1212
P3 = 1.6151
P4 = 1.4066
P5 = -0.3040
P6 = 0.4843
P7 = -0.3285
P8 = -2.0749
A = P1 * (P_crit ** P2) * (G_Mlbs ** (P5 + (P7 * P_crit)))
C = P3 * (P_crit ** P4) * (G_Mlbs ** (P6 + (P8 * P_crit)))
q_crit_heat_flux_MBtu = (A - xe_array[0]) / (C + ((xe_array[1:] - xe_array[0]) / q_heat_flux_MBtu))
q_crit_heat_flux = q_crit_heat_flux_MBtu * (1 / 3.41) * 1000000 * 10.7639
DNBR = q_crit_heat_flux / q_heat_flux_DNBR
plt.figure(4)
plt.plot(DNBR, z_arrayplots_DNBR)
plt.xlabel("Departure from Nucleate Boiling Ratio")
plt.ylabel("Height z [m]")
plt.savefig("DNBR.png", dpi=600)
plt.figure(5)
plt.plot(P_array, z_arrayplots)
plt.xlabel('Pressure [MPa]')
plt.ylabel('Height z [m]')
plt.savefig("Pressure.png", dpi=600)
plt.figure(6)
plt.plot(T_f_array, z_arrayplots)
plt.xlabel('Temperature [K]')
plt.ylabel('Height z [m]')
plt.savefig("TempBulk.png", dpi=600)
plt.figure(7)
plt.plot(T_F_Outer, z_arrayplots, label='Fuel Outer Surface Temp')
plt.plot(T_F_Center, z_arrayplots, label='Fuel Centerline Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempFuel.png", dpi=600)
tempdifference = T_C_Outer - T_f_array
print("Max clad vs bulk difference is " + str(np.amax(tempdifference)) + " K")
print("Max coolant temp is " + str(np.amax(T_f_array)) + " K")
print("Min coolant temp is " + str(np.amin(T_f_array)) + " K")
print("Max clad inner temp is " + str(np.amax(T_C_Inner)) + " K")
print("Max clad outer temp is " + str(np.amax(T_C_Outer)) + " K")
print("min clad outer temp is " + str(np.amin(T_C_Outer)) + " K")
print("Max fuel temp is " + str(np.amax(T_F_Center)) + " K")
print("Max fuel outer temp is " + str(np.amax(T_F_Outer)) + " K")
print("Min fuel outer temp is " + str(np.amin(T_F_Outer)) + " K")
print("Max centerline temp occurs at z = " + str(z_CL_max) + "m")
print("Max clad temp occurs at z = " + str(z_Clad_max) + "m")
MDNBR = np.amin(DNBR)
print("MDNBR is " + str(MDNBR))
plt.show()
if case == 2: # BWR
H = 3.8 # m
He = 3.8 # m
Pitch = 1.63 * 10 ** (-2) # m
Gap_t = 0.0001 # m
D_fuel = 0.0104 # m
k_gap = 0.25 # W/m-K
k_c = 21.5 # W/m-K
k_fuel = 3.6 # W/m-K
T0 = float(274 + 273.15) # K
q0_prime = float(410 * 10 ** (2)) # W/m
P0 = 7.5 # MPa
MF = float(2290) # kg/m^2-s
D_rod = .0123 # m
R_rod = D_rod / 2
R_fuel = D_fuel / 2
R_gap = R_fuel + Gap_t
R_clad = R_rod
Clad_t = D_rod - D_fuel - Gap_t # m
h0_enthalpy = (IAPWS97(T=T0, P=P0).h) * 10 ** (3)
T_sat0 = IAPWS97(P=P0, x=0).T
g = 9.81 # m/s
# geometry properties
heated_p = math.pi * D_rod
wetted_p = math.pi * D_rod
A_f = (Pitch ** 2) - ((1 / 4) * math.pi * (D_rod ** 2))
D_H = (4 * A_f) / heated_p
# grid setup
grid_points = 100
dz = H / grid_points
z_array = np.arange(0, H, dz)
z_arrayplots = np.arange(-H / 2, H / 2, dz)
q_HeatFluxList = []
# defining array of q'' values in list
for z in z_array:
heat_fluxA = (q0_prime / (math.pi * D_rod)) * math.sin(math.pi * (z / He))
q_HeatFluxList.append(heat_fluxA)
q_heat_flux = np.array(q_HeatFluxList)
q_prime = np.zeros(len(z_array))
for i in range(0, len(z_array)):
q_prime[i] = q0_prime * math.sin((np.pi * z_array[i]) / He)
# defining array of h values
h_enthalpy_list = []
h_enthalpy_prefactor = ((heated_p * q0_prime * H) / (A_f * MF * (math.pi ** 2) * D_rod))
for z in z_array:
h_enthalpy = (-h_enthalpy_prefactor * math.cos(math.pi * (z / He))) + h_enthalpy_prefactor + h0_enthalpy
h_enthalpy_list.append(h_enthalpy)
h_enthalpy_array_J = np.array(h_enthalpy_list)
h_enthalpy_array = h_enthalpy_array_J * 10 ** (-3)
P_array = np.zeros(len(z_array))
P_array[0] = P0
T_sat = np.zeros(len(z_array))
T_sat[0] = T_sat0
T_f_array = np.zeros(len(z_array))
T_f_array[0] = T0
Re = np.zeros(len(z_array))
Re_f = np.zeros(len(z_array))
Pr = np.zeros(len(z_array))
k_fluid = np.zeros(len(z_array))
x_array = np.zeros(len(z_array))
xe_array = np.zeros(len(z_array))
dxe_array = np.zeros(len(z_array))
frictional = np.zeros(len(z_array))
gravitational = np.zeros(len(z_array))
compressibility = np.zeros(len(z_array))
alpha_array = np.zeros(len(z_array))
# Pressure Loop BWR
dp = 0.001
for i in range(0, len(z_array) - 1):
rho_f = IAPWS97(P=P_array[i], x=0).rho
rho_m = IAPWS97(P=P_array[i], x=xe_array[i]).rho
vf = IAPWS97(P=P_array[i], x=0).v
vg = IAPWS97(P=P_array[i], x=1).v
vfg = vg - vf
hf_enthalpy = IAPWS97(P=P_array[i], x=0).h
hg_enthalpy = IAPWS97(P=P_array[i], x=1).h
hfg = hg_enthalpy - hf_enthalpy
hfg_sat = IAPWS97(P=P0, x=1).h - IAPWS97(P=P0, x=0).h
# vf = IAPWS97(P=P_array[i], x=0).v
# vg_sat = IAPWS97(P=P_array[i], x=1).v
hf_in = IAPWS97(P=P0, T=T0).h
muf = (IAPWS97(P=P_array[i], x=0).mu) * 10 ** (-6)
mum = (IAPWS97(P=P_array[i], x=xe_array[i]).mu) * 10 ** (-6)
mug = (IAPWS97(P=P_array[i], x=1).mu) * 10 ** (-6)
k_fluid[i] = IAPWS97(P=P_array[i], T=T_f_array[i]).k
Pr[i] = IAPWS97(P=P_array[i], h=h_enthalpy_array[i]).Liquid.Prandt
xe_in = (hf_in - hf_enthalpy) / (hfg)
vf_sat = IAPWS97(P=P_array[i], x=xe_array[i])
# vapor quality
if xe_array[i] <= 0: # single phase
Re1p = MF * D_rod / muf
f1p = 0.316 * Re1p ** (-.25)
dp = -(.5 * f1p * MF * 2 * heated_p / (rho_f * A_f) + g / rho_f) * dz
x_array[i] = 0
dxe_array[i] = q0_prime * np.sin(np.pi * z_array[i] / H) / (MF * A_f * hfg_sat) * dz
xe_array = xe_array[i - 1] + dxe_array[i]
# P_array[i]=P[i-1]+dp1p
elif xe_array[i] > 0 and xe_array[i] < 1: # 2 phase
Re2p = MF * D_rod / mum
f2p = 0.046 * Re2p ** (-.2) * (muf / mum ** (-.2))
dp2p = (-MF ** 2 * vfg * dxe_array[i] + .5 * f2p * MF ** 2 * heated_p / (rho_m * A_f) + g * rho_m) * dz
xe_array[i] = (h_enthalpy_array[i] - hf_enthalpy) / (hg_enthalpy - hf_enthalpy)
# Void Fraction
if xe_array[i] <= 0:
alpha_array[0]
elif xe_array[i] > 0 and xe_array[i] < 1:
x_array[i] = xe_array[i]
vfg_sat = vg - vf
rho_m = (vf_sat + vfg_sat * x_array) ** (-1)
rhof = 1 / vf
rhog = 1 / vg
void = (rho_m - rhof) / (rhog - rhof)
alpha_array[void]
print("Void fraction is " + str(np.amax(alpha_array)))
if xe_array[i] <= 0:
alpha_array[0]
elif xe_array[i] > 0 and xe_array[i] < 1:
x_array[i] = xe_array[i]
vfg_sat = vg - vf
rho_m = (vf_sat + vfg_sat * x_array) ** (-1)
rhof = 1 / vf
rhog = 1 / vg
void = (rho_m - rhof) / (rhog - rhof)
alpha_array[void]
print("Void fraction is " + str(np.amax(alpha_array)))
rho_m = 1 / ((x_array[i] * vg) + ((1 - x_array[i]) * vf))
mu_m = 1 / ((x_array[i] / mug) + ((1 - x_array[i]) / muf))
Re[i] = (MF * D_H) / (mu_m * 10 ** 6) # convert mu to Pa/s
f = 0.079 * (Re[i] ** -0.25) * (mu_m / muf)
Tau = (1 / 2) * f * ((MF ** 2) / rho_m)
Re_f[i] = Re[i]
vf_plus_dP = IAPWS97(P=P_array[i] + dp, x=xe_array[i]).v
vf_minus_dP = IAPWS97(P=P_array[i] - dp, x=xe_array[i]).v
ddP_vf = (vf_plus_dP - vf_minus_dP) / (2 * (dp * 10 ** 6))
frictional[i] = (Tau * wetted_p) / A_f
gravitational[i] = g * rho_f
compressibility[i] = (MF ** 2) * (ddP_vf)
dPdz_num = (frictional[i] + gravitational[i]) # Pa/m
dPdz_denom = 1 + compressibility[i] # Pa/m
dPdz = -dPdz_num / dPdz_denom # Pa/m
P_array[i + 1] = P_array[i] + ((dPdz * dz) * 10 ** (-6))
T_f_array[i + 1] = IAPWS97(P=P_array[i + 1], h=h_enthalpy_array[i + 1]).T
T_sat[i + 1] = IAPWS97(P=P_array[i + 1], x=0).T
# final calc for final value of quality and void fraction because loop stops before these
hf_final = IAPWS97(P=P_array[-1], x=0).h
hg_final = IAPWS97(P=P_array[-1], x=1).h
muf_final = (IAPWS97(P=P_array[-1], x=0).mu) * 10 ** (-6)
mug_final = (IAPWS97(P=P_array[-1], x=1).mu) * 10 ** (-6)
k_fluid[-1] = IAPWS97(P=P_array[-1], T=T_f_array[-1]).k
xe_array[-1] = (h_enthalpy_array[-1] - hf_final) / (hg_final - hf_final)
# fuel and clad temps
T_C_Outer = np.zeros(len(z_array))
mu_m_final = 1 / ((x_array[-1] / mug_final) + ((1 - x_array[-1]) / muf_final))
Re_f[-1] = (MF * D_H) / (muf_final * 10 ** 6)
Pr[-1] = IAPWS97(P=P_array[i], h=h_enthalpy_array[i]).Liquid.Prandt
h_HT = 0.023 * (Re_f[0] ** 0.8) * (Pr[0] ** 0.4) * (k_fluid[0] / D_H)
T_C_Outer[0] = (q_heat_flux[0] + (h_HT * T_f_array[0])) / h_HT
for i in range(0, len(z_array) - 1):
h_HT = 0.023 * (Re_f[i + 1] ** 0.8) * (Pr[i + 1] ** 0.4) * (k_fluid[i + 1] / D_H)
T_C_Outer[i + 1] = (q_heat_flux[i + 1] + (h_HT * T_f_array[i + 1])) / h_HT
q_triple_prime = (q_prime * 4) / (np.pi * (D_fuel ** 2))
T_C_Inner = np.zeros(len(z_array))
T_F_Outer = np.zeros(len(z_array))
T_F_Center = np.zeros(len(z_array))
for i in range(0, len(z_array)):
C1 = -((q0_prime * R_clad) / (k_c * heated_p)) * np.sin(np.pi * (z_array[i] / H))
C2 = T_C_Outer[i] - (C1 * np.log(R_clad))
T_C_Inner[i] = (C1 * np.log(R_gap)) + C2
C3 = (k_c / k_gap) * C1
C4 = T_C_Inner[i] - (C3 * np.log(R_gap))
T_F_Outer[i] = (C3 * np.log(R_fuel)) + C4
C6 = T_F_Outer[i] + ((q_triple_prime[i] * (R_fuel ** 2)) / (4 * k_fuel))
T_F_Center[i] = C6
CL_max = np.amax(T_F_Center)
index = np.where(T_F_Center == CL_max)
z_CL_max = z_array[index]
plt.figure(1)
plt.plot(T_C_Outer, z_arrayplots, label='Clad Outer Surface Temp')
plt.plot(T_C_Inner, z_arrayplots, label='Clad Inner Surface Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempCladBWR.png", dpi=600)
plt.figure(2)
plt.plot(T_C_Outer, z_arrayplots, label='Clad Outer Surface Temp')
plt.plot(T_C_Inner, z_arrayplots, label='Clad Inner Surface Temp')
plt.plot(T_F_Outer, z_arrayplots, label='Fuel Outer Surface Temp')
plt.plot(T_F_Center, z_arrayplots, label='Fuel Centerline Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempFuelAndCladBWR.png", dpi=600)
# radial calcs
T_array_A = [T_F_Center[25], T_F_Outer[25], T_C_Inner[25], T_C_Outer[25]]
T_array_B = [T_F_Center[49], T_F_Outer[49], T_C_Inner[49], T_C_Outer[49]]
T_array_C = [T_F_Center[53], T_F_Outer[53], T_C_Inner[53], T_C_Outer[53]]
r_array = [0, R_fuel, R_gap, R_clad]
plt.figure(3)
plt.plot(r_array, T_array_A, label='z = -H/4 = -0.9 m')
plt.plot(r_array, T_array_B, label='z = 0 m')
plt.plot(r_array, T_array_C, '--', label='z = zmax = 0.108 m')
plt.legend(loc='upper left')
plt.ylabel("Temperature [K]")
plt.xlabel("Radius r [m]")
plt.savefig("TempRadialBWR.png", dpi=600)
# critical heat flux and DNBR
P_array_DNBR = np.delete(P_array, 0)
q_heat_flux_DNBR = np.delete(q_heat_flux, 0)
z_arrayplots_DNBR = np.delete(z_arrayplots, 0)
G_Mlbs = MF * (((2.20462 * 10 ** (-6)) * 3600) / 10.7639)
q_heat_flux_MBtu = q_heat_flux[1:] * 3.41 * (1 / 1000000) * (1 / 10.7639)
P_c = 22.064 # https://nuclearstreet.com/nuclear-power-plants/w/nuclear_power_plants/features-of-pressurized-water-reactors
P_crit = P_array_DNBR / P_c
P1 = 0.5328
P2 = 0.1212
P3 = 1.6151
P4 = 1.4066
P5 = -0.3040
P6 = 0.4843
P7 = -0.3285
P8 = -2.0749
A = P1 * (P_crit ** P2) * (G_Mlbs ** (P5 + (P7 * P_crit)))
C = P3 * (P_crit ** P4) * (G_Mlbs ** (P6 + (P8 * P_crit)))
q_crit_heat_flux_MBtu = (A - xe_array[0]) / (C + ((xe_array[1:] - xe_array[0]) / q_heat_flux_MBtu))
q_crit_heat_flux = q_crit_heat_flux_MBtu * (1 / 3.41) * 1000000 * 10.7639
DNBR = q_crit_heat_flux / q_heat_flux_DNBR
plt.figure(4)
plt.plot(DNBR, z_arrayplots_DNBR)
plt.xlabel("Onset of Nucleate Boiling Ratio")
plt.ylabel("Height z [m]")
plt.title("Onset of Nucleate Boiling Ratio versus Height")
plt.savefig("ONBR.png", dpi=600)
plt.figure(5)
plt.plot(P_array, z_arrayplots)
plt.xlabel('Pressure [MPa]')
plt.ylabel('Height z [m]')
plt.title('Pressure versus Height')
plt.savefig("PressureBWR.png", dpi=600)
plt.figure(6)
plt.plot(T_f_array, z_arrayplots)
plt.xlabel('Temperature [K]')
plt.ylabel('Height z [m]')
plt.title('Coolant Temperature vs Height')
plt.savefig("TempBulkBWR.png", dpi=600)
plt.figure(7)
plt.plot(T_F_Outer, z_arrayplots, label='Fuel Outer Surface Temp')
plt.plot(T_F_Center, z_arrayplots, label='Fuel Centerline Temp')
plt.legend(loc='upper left')
plt.xlabel("Temperature [K]")
plt.ylabel("Height z [m]")
plt.savefig("TempFuelBWR.png", dpi=600)
# density
plt.figure(8)
plt.plot(Density, z_arrayplots, label='Density')
plt.legend(loc='upper left')
plt.xlabel("Pressure [mPa]")
plt.ylabel("Height z [m]")
plt.savefig("Density", dpi=600)
# quality
plt.figure(9)
plt.plot(x, z_arrayplots, label='Quality')
plt.plot(xe, z_arrayplots, label='Quality')
plt.legend(loc='upper left')
plt.xlabel("Quality")
plt.ylabel("Height z [m]")
plt.savefig("Quality", dpi=600)
# void
plt.figure(10)
plt.plot(alpha, z_arrayplots, label='Void Fraction')
plt.legend(loc='upper left')
plt.xlabel("Void Fraction")
plt.ylabel("Height z [m]")
plt.savefig("Void Fraction", dpi=600)
tempdifference = T_C_Outer - T_f_array
print("Max clad vs bulk difference is " + str(np.amax(tempdifference)) + " C")
print("Max coolant temp is " + str(np.amax(T_f_array) - 273.15) + " C")
print("Max coolant temp is " + str(np.amax(T_f_array)) + " K")
print("Max clad temp is " + str(np.amax(T_C_Inner) - 273.15) + " C")
print("Max clad temp is " + str(np.amax(T_C_Inner)) + " K")
print("Max fuel temp is " + str(np.amax(T_F_Center) - 273.15) + " C")
print("Max fuel temp is " + str(np.amax(T_F_Center)) + " K")
print("Max fuel temp is " + str(np.amax(T_F_Outer) - 273.15) + " C")
print("Max fuel temp is " + str(np.amax(T_F_Outer)) + " K")
print("Max centerline temp occurs at z = " + str(z_CL_max) + "m")
MDNBR = np.amin(DNBR)
print("MDNBR is " + str(MDNBR))
PWR 和 BWR 情况的区别在于，BWR 有两相水（液体和蒸汽），而 PWR 只有一相水（液体）。
IAPWS97 库要求用户在有多个阶段时指定他们想要的阶段：
IAPWS97(T=T0, P=P0).h # if only one phase
应该是
IAPWS97(T=T0, P=P0).Liquid.h
或
IAPWS97(T=T0, P=P0).Vapor.h