crp/MSRE
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

add dynamic model

Browse source
This commit is contained in:
lukelabrie 2023-08-07 16:52:55 +02:00
parent 67ca89f66b
commit 5930aeceec
3 changed files with 577 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

249
dynamic_model/main.py Normal file
View file

@ -0,0 +1,249 @@
from parameters import *
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import copy
def dydtMSRE(t,y,delays,rho_ext):
'''
Returns derivative of state vector y; y' or dy/dt, of the MSRE system.
The y vector contains the following:
T_in_rc = Inlet temperature (°C) of radiator coolant, will be equal to
the outlet temperature of the heat exchanger (T_hc4) plus
relevant time delay
T_out_rc = Outlet temperature (°C) of radiator coolant
T_in_air = Inlet temperature (°C) of air in radiator
T_out_air = Outlet temperature (°C) of air in radiator
T_in_hf = Inlet temperature (°C) of heat exchanger fuel, will be equal
to the outlet temperature of the core (T_cf2) plus relevant
time delay
T_hf* = Temperature (°C) of heat exchanger fuel node *
T_ht* = Temperature (°C) of heat exchanger tube node *
T_hc* = Temperature (°C) of heat exchanger coolant node *
T_in_cf = Inlet temperature (°C) of core fuel, will be equal to the
outlet temperature of the heat exchanger (T_hf4) plus
relevant time dela
S = neutro source perturbation term
rho_fb = feedback reactivity (from fuel and graphite temperatures)
rho_ext = external reactivity (reactivity insertion)
rho_tot = total reactivity = rho_0 + rho_fb + rho_ext (rho_0 =
steady-state reactivity, constant)
n = neutron density n(t)
C* = precursor concentration of group *
T_cg = Temperature (°C) of core graphite node
T_cf* = Temperature (°C) of core fuel node *
*_delay = Parameter at time t = t - delay
Other parameters are defined in parameters.py
'''
# unpack state variables
T_out_rc, T_out_air, T_hf1, T_hf2, T_hf3, T_hf4, T_ht1, T_ht2, T_hc1, \
T_hc2, T_hc3, T_hc4, n, C1, C2, C3, C4, C5, C6, T_cg, T_cf1, T_cf2 = y
# delay terms
T_out_rc_delay, T_hc4_delay, T_cf2_delay, T_hf4_delay, C1_delay, \
C2_delay, C3_delay, C4_delay, C5_delay, C6_delay = delays
# reactivity
rho = (a_f/2)*((-T0_f1+T_cf1)+(-T0_f2+T_cf2)) + a_g*(-T0_g1+T_cg) + rho_ext
# derivatives
dydt = [
(W_rp/mn_rp)*(T_hc4_delay-T_out_rc) + (hA_rpn/mcp_rpn)*(T_out_air-T_out_rc), # T_out_rc # T_out_rc
-((W_rs/mn_rs)+(hA_rsn/mcp_rsn))*T_out_air + (hA_rsn/mcp_rsn)*T_out_rc + (W_rs/mn_rs)*Trs_in, # T_out_air
-((W_p/mn_p)+(hA_pn/mcp_pn))*T_hf1 + (hA_pn/mcp_pn)*T_ht1 + (W_p/mn_p)*T_cf2_delay, # T_hf1
(W_p/mn_p)*(T_hf1-T_hf2) + (hA_pn/mcp_pn)*(T_ht1-T_hf1), # T_hf2
-((W_p/mn_p)+(hA_pn/mcp_pn))*T_hf3 + (hA_pn/mcp_pn)*T_ht2 + (W_p/mn_p)*T_hf2, # T_hf3
(W_p/mn_p)*(T_hf3-T_hf4) + (hA_pn/mcp_pn)*(T_ht2-T_hf3), # T_hf4
(2*hA_pn/mcp_tn)*(T_hf1-T_ht1) + (2*hA_sn/mcp_tn)*(T_hc3-T_ht1), # T_ht1
(2*hA_pn/mcp_tn)*(T_hf3-T_ht2) + (2*hA_sn/mcp_tn)*(T_hc1-T_ht2), # T_ht2
-((W_s/mn_s)+(hA_sn/mcp_sn))*T_hc1 + (hA_sn/mcp_sn)*T_ht2 + (W_s/mn_s)*T_out_rc_delay, # T_hc1
(W_s/mn_s)*(T_hc1-T_hc2) + (hA_sn/mcp_sn)*(T_ht2-T_hc1), # T_hc2
-((W_s/mn_s)+(hA_sn/mcp_sn))*T_hc3 + (hA_sn/mcp_sn)*T_ht1 + (W_s/mn_s)*T_hc2, # T_hc3
(W_s/mn_s)*(T_hc3-T_hc4) + (hA_sn/mcp_sn)*(T_ht1-T_hc3), # T_hc4
(rho-beta_t)*n/Lam+lam[0]*C1+lam[1]*C2+lam[2]*C3+lam[3]*C4+lam[4]*C5+lam[5]*C6, # n (no source insertion)
n*beta[0]/Lam-lam[0]*C1-C1/tau_c+C1_delay*np.exp(-lam[0]*tau_l)/tau_c, # C1
n*beta[1]/Lam-lam[1]*C2-C2/tau_c+C2_delay*np.exp(-lam[1]*tau_l)/tau_c, # C2
n*beta[2]/Lam-lam[2]*C3-C3/tau_c+C3_delay*np.exp(-lam[2]*tau_l)/tau_c, # C3
n*beta[3]/Lam-lam[3]*C4-C4/tau_c+C4_delay*np.exp(-lam[3]*tau_l)/tau_c, # C4
n*beta[4]/Lam-lam[4]*C5-C5/tau_c+C5_delay*np.exp(-lam[4]*tau_l)/tau_c, # C5
n*beta[5]/Lam-lam[5]*C6-C6/tau_c+C6_delay*np.exp(-lam[5]*tau_l)/tau_c, # C6
(hA_fg/mcp_g1)*(T_cf1 - T_cg) + k_g*P*n/mcp_g1, # T_cg
W_f/mn_f*(T_hf4_delay-T_cf1) + (k_f1*P*n/mcp_f1) + (hA_fg*k_1*(T_cg - T_cf1)/mcp_f1), # T_cf1
W_f/mn_f*(T_cf1 - T_cf2) + (k_f2*P*n/mcp_f2) + (hA_fg*k_2*(T_cg - T_cf1)/mcp_f2) # T_cf2
]
return dydt
def get_tIdx(t,tao,timeVec):
'''
Returns index of time t = t-tau
'''
td = t-tao
diff_min = 999999.9999999
idx = 0
for t in enumerate(timeVec):
diff = abs(td-t[1])
if (diff<diff_min):
diff_min = abs(td-t[1])
idx = t[0]
return idx, timeVec[idx]-td
def main():
'''
Sets initial conditions and calls the solver
'''
# initial conditions
y0 = [T0_rp, T0_rs, T0_p1,T0_p2, T0_p3, T0_p4, T0_t1, T0_t2, T0_s1, T0_s2,
T0_s3, T0_s4, n_frac0, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5],
T0_g1, T0_f1, T0_f2]
# initial delay terms
d0 = [T0_rp, T0_rs, T0_f2, T0_p4, C0[0], C0[1], C0[2], C0[3], C0[4], C0[5]]
# solver
backend = 'dopri5'
r = ode(dydtMSRE).set_integrator(backend,max_step=0.10)
sol_interim = []
def solout(t, y):
sol_interim.append([t, *y])
r.set_solout(solout)
# timing parameters
t0 = 0.0
t_start = t0
t_stop = 5000
# solution
sol = []
# step-reactivity insertion
t_insert = 2500.00
insert = 1.0e-4
if (t_start>=t_insert):
rho_ext = insert
else:
rho_ext = 0.0
# delay parameters
d_terms = []
derivs = [dydtMSRE(t0,y0,d0,rho_ext)]
# takes one step at a time and then accounts for delay terms
i = 0
while (t_start < t_stop):
# take one step
if (i == 0):
t_start = t0
r.set_initial_value(y0,t0).set_f_params(d0,0.0)
r.integrate(1.0)
sol.append(sol_interim[0])
sol.append(sol_interim[1])
else:
t_start = sol[-1][0]
if (t_start>=t_insert):
rho_ext = insert
else:
rho_ext = 0.0
r.set_initial_value(y_next,t_start).set_f_params(d_new,rho_ext)
r.integrate(t_start+1.0)
sol.append(sol_interim[1])
derivs.append(dydtMSRE(t_start,y_next,d_new,rho_ext))
# account for delays, linear interpolation for time differences
# core fuel inlet
d_new = [0]*10
idx_cf_in = 0
dt_cf = 0.0
if (t_start > tau_hx_c):
idx_cf_in, dt_cf = get_tIdx(t_start,tau_hx_c,[s[0] for s in sol])
d_new[3] = sol[idx_cf_in][6] + dt_cf*derivs[idx_cf_in][5]
# heat exchanger fuel inlet
idx_hf_in = 0
dt_hf = 0.0
if (t_start > tau_c_hx):
idx_hf_in, dt_hf = get_tIdx(t_start,tau_c_hx,[s[0] for s in sol])
d_new[2] = sol[idx_hf_in][22] + dt_hf*derivs[idx_hf_in][21]
# heat exchanger coolant inlet
idx_hc_in = 0
dt_hc = 0.0
if (t_start > tau_r_hx):
idx_hc_in, dt_hc = get_tIdx(t_start,tau_r_hx,[s[0] for s in sol])
d_new[0] = sol[idx_hc_in][1] + dt_hc*derivs[idx_hc_in][0]
# radiator coolant inlet
idx_rc_in = 0
dt_rc = 0.0
if (t_start > tau_hx_r):
idx_rc_in, dt_rc = get_tIdx(t_start,tau_hx_r,[s[0] for s in sol])
d_new[1] = sol[idx_rc_in][12] + dt_rc*derivs[idx_rc_in][11]
# precursors
idx_c = 0
dt_c = 0.0
if (t_start > tau_l):
idx_c, dt_c = get_tIdx(t_start,tau_l,[s[0] for s in sol])
d_new[4] = sol[idx_c][14] + dt_c*derivs[idx_c][13]
d_new[5] = sol[idx_c][15] + dt_c*derivs[idx_c][14]
d_new[6] = sol[idx_c][16] + dt_c*derivs[idx_c][15]
d_new[7] = sol[idx_c][17] + dt_c*derivs[idx_c][16]
d_new[8] = sol[idx_c][18] + dt_c*derivs[idx_c][17]
d_new[9] = sol[idx_c][19] + dt_c*derivs[idx_c][18]
d_terms.append(d_new)
# initial condiiton for next step
y_next = sol[-1][1:]
# empty interim solution
sol_interim = []
# display progress
print(f"{t_start}")
i += 1
# plot single parameter
of_interest = 14
plt.plot([s[0] for s in sol],[s[of_interest] for s in sol])
plt.show()
# check delay behavior
# for i in range(len(sol)-1):
# print(f"t: {sol[i][0]}, hf4: {sol[i][6]}, c1_delay: {d_terms[i][3]}")
# write output data
output_filename = f"sim_out_{t_stop}_{P}"
results = open(output_filename,'w+')
for k in range(len(sol)):
for col in range(len(sol[0])):
results.write(f"{sol[k][col]} ")
results.write("\n")
return None
main()

234
dynamic_model/parameters.py Normal file
View file

@ -0,0 +1,234 @@
import numpy as np
import pandas as pd
import math
pi = math.pi
# Perturbations
# SOURCE INSERTION
# No source insertion
sourcedata = np.array([0, 0, 0])
sourcetime = np.array([0, 50, 100])
# % 1 (n/no)/s for 10 seconds
# sourcedata = np.array([0, 10, 0])
# sourcetime = np.array([0, 10, 20])
source = pd.Series(sourcedata, index=sourcetime)
# REACTIVITY INSERTION
# No reactivity insertion
simtime = 10
reactdata = np.array([0, 5E-4])
reacttime = np.array([0, 2500])
# Periodic 60 PCM for 50 seconds
# simtime = 500
# periodic = np.array([[0, 0], [50, 6e-4], [100, 0], [150, -6e-4], [200, 0], [250, 6e-4], [300, 0], [350, -6e-4], [400, 0]])
# reactdata = periodic[:, 1]
# reacttime = periodic[:, 0]
# Step up 60 pcm
# simtime = 1000
# reactdata = np.array([0, 6e-3])
# reacttime = np.array([0, 300])
# # Step down -60 pcm for 10 sec
# simtime = 100
# reactdata = np.array([0, -6e-4])
# reacttime = np.array([0, 50])
# # Pulse 600 pcm for 0.1 sec
# simtime = 30
# reactdata = np.array([0, 6e-3, 0])
# reacttime = np.array([0, 10, 10.1])
react = pd.Series(reactdata, index=reacttime)
ts_max = 1e-1 # maximum timestep (s)
# NEUTRONICS DATA
tau_l = 16.73 # ORNL-TM-0728 %16.44; % (s)
tau_c = 8.46 # ORNL-TM-0728 %8.460; % (s)
P = 8 # Thermal Power in MW ORNL-TM-1070, p.2
n_frac0 = 5 # initial fractional neutron density n/n0 (n/cm^3/s)
Lam = 2.400E-04 # mean generation time ORNL-TM-1070 p.15 U235
# Lam = 4.0E-04; # mean generation time ORNL-TM-1070 p.15 U233
lam = np.array([1.240E-02, 3.05E-02, 1.11E-01, 3.01E-01, 1.140E+00, 3.014E+00])
beta = np.array([0.000223, 0.001457, 0.001307, 0.002628, 0.000766, 0.00023]) # U235
# beta = np.array([0.00023, 0.00079, 0.00067, 0.00073, 0.00013, 0.00009]) # U233
beta_t = np.sum(beta) # total delayed neutron fraction MSRE
rho_0 = beta_t-sum(np.divide(beta,1+np.divide(1-np.exp(-lam*tau_l),lam*tau_c))) # reactivity change in going from stationary to circulating fuel
C0 = beta / Lam * (1.0 / (lam - (np.exp(-lam * tau_l) - 1.0) / tau_c))
# Feedback co-efficients
a_f = -8.71E-05 # U235 (drho/°C) fuel salt temperature-reactivity feedback coefficient ORNL-TM-1647 p.3 % -5.904E-05; % ORNL-TM-0728 p. 101 %
a_g = -6.66E-05 # U235 (drho/°C) graphite temperature-reactivity feedback coefficient ORNL-TM-1647 p.3 % -6.624E-05; % ORNL-TM-0728 p.101
# CORE HEAT TRANSFER PARAMETERS
# FUEL PARAMETERS - DONE
vdot_f = 7.5708E-02 # ORNL-TM-0728 % 7.571e-2; % vol. flow rate (m^3/s) ORNL-TM-1647 p.3, ORNL-TM-0728 p.12
rho_f = 2.14647E+03 # (partially enriched U-235)ORNL-TM-0728 p.8 2.243E+03; % (Th-U) density of fuel salt (kg/m^3) ORNL-TM-0728 p.8
W_f = 1.623879934566580e+02 # 1.83085e+02;%vdot_f*rho_f; % 182.78; % calcd from m_dot*cp*delT=P; vdot_f*rho_f; % fuel flow rate (kg/s)
# tau_f_c = tau_c; % ORNL-TM-0728 % 8.45; % transit time of fuel in core (s) ORNL-TM-1070 p.15, TDAMSRE p.5
m_f = W_f * tau_c # fuel mass in core (kg)
nn_f = 2 # number of fuel nodes in core model
mn_f = m_f / nn_f # fuel mass per node (kg)
# cp_f = 4.2*9/5; % (MJ/deg-C) total fuel heat capacity TDAMSRE p.5
scp_f = 1.9665E-3 # specific heat capacity of fuel salt (MJ/kg-C) ORNL-TM-0728 p.8
# Core Upflow - DONE
v_g = 1.95386 # graphite volume(m^3) ORNL-TM-0728 p. 101
rho_g = 1.860E3 # graphite density (kg/m^3) ORNL-3812 p.77, ORNL-TM-0728 p.87
m_g = v_g * rho_g # graphite mass (kg)
cp_g = 3.6 * 9 / 5 # TDAMSRE p.5 graphite total heat capacity (MW-s/C) ORNL-TM-1647 p.3
scp_g = 1.773E-3 # cp_g/m_g; % graphite specific heat capacity (MW-s/kg-C) ORNL-TM-1647 p.3
mcp_g1 = m_g * scp_g # (mass of material x heat capacity of material) of graphite per lump (MW-s/°C)
mcp_f1 = mn_f * scp_f # (mass of material x heat capacity of material) of fuel salt per lump (MW-s/°C)
mcp_f2 = mn_f * scp_f # (mass of material x heat capacity of material) of fuel salt per lump (MW-s/°C)
hA_fg = 0.02 * 9 / 5 # (fuel to graphite heat transfer coeff x heat transfer area) (MW/°C) ORNL-TM-1647 p.3, TDAMSRE p.5
k_g = 0.07 # fraction of total power generated in the graphite ORNL-TM-0728 p.9
k_1 = 0.5 # fraction of heat transferred from graphite which goes to the first fuel lump
k_2 = 0.5 # fraction of heat transferred from graphite which goes to the second fuel lump
k_f = 0.93 # fraction of heat generated in fuel - that generated in the external loop ORNL-TM-0728 p.9
k_f1 = k_f / nn_f # fraction of total power generated in lump f1
k_f2 = k_f / nn_f # fraction of total power generated in lump f2
# New node for power deposited in fuel outside the core
k_out = 1 - (k_g + k_f) # fraction of power generated in fuel in external loop ORNL-TM-0728 p.9
m_out = W_f # (kg) Mass of node such that resident time is 1 sec (W_f needs to be defined)
# Initial conditions - DONE
Tf_in = 6.3222E+02 # in °C ORNL-TM-1647 p.2
T0_f2 = 6.5727E+02 # 6.5444E+02; % in °C 6.461904761904777e+02; ORNL-TM-1647 p.2
T0_f1 = Tf_in + (T0_f2 - Tf_in) / 2 # 6.405952380952389e+02; in °C
T0_g1 = T0_f1 + (k_g * P / hA_fg) # 6.589285714285924e+02; in °C
# T0_out = k_out * P / m_out / scp_f + T0_f2 # in °C (scp_f needs to be defined)
# Heat Exchanger - DONE
# Geometry
d_he = 16 # (in) he diameter ORNL-TM-0728 p. 164
h_he = 72 # (in) active height % 96; %(in) he height ORNL-TM-0728 p. 164
od_tube = 0.5 # (in) coolant tube OD ORNL-TM-0728 p. 164
id_tube = od_tube - 2 * 0.042 # (in) coolant tube ID ORNL-TM-0728 p. 164
n_tube = 159 # number of coolant tubes ORNL-TM-0728 p. 164
a_tube = 254 * 144 # (in^2) total area of tubes ORNL-TM-0728 p. 164
l_tube = a_tube / n_tube / (np.pi * od_tube) # (in) tube length
v_tube = n_tube * np.pi * (od_tube / 2) ** 2 * l_tube # (in^3) hx shell volume occupied by tubes
v_cool = n_tube * np.pi * (id_tube / 2) ** 2 * l_tube # (in^3) hx volume occupied by coolant
v_he = (d_he / 2) ** 2 * np.pi * h_he # (in^3) volume of heat exchanger shell
v_he_fuel = v_he - v_tube # (in^3) volume available to fuel in shell
# Unit conversions
in_m = 1.63871e-5 # 1 cubic inch = 1.63871e-5 cubic meters
# PRIMARY FLOW PARAMETERS - DONE
W_p = W_f # fuel flow rate (kg/s)
m_p = v_he_fuel * in_m * rho_f # fuel mass in PHE (kg)
nn_p = 4 # number of fuel nodes in PHE
mn_p = m_p / nn_p # fuel mass per node (kg)
cp_p = scp_f # fuel heat capacity (MJ/(kg-C))
# SECONDARY FLOW PARAMETERS - DONE
vdot_s = 5.36265E-02 # ORNL-TM-0728 p. 164 % 5.236E-02; % coolant volume flow rate (m^3/s) ORNL-TM-1647 p.3
rho_s = 1.922e3 # coolant salt density (kg/m^3) ORNL-TM-0728 p.8
W_s = 1.005793369810108e+02 # vdot_s*rho_s; % calcd from mdot*cp*delT; vdot_s*rho_s; % coolant flow rate (kg/s) ORNL-TM-1647 p.3
m_s = v_cool * in_m * rho_s # coolant mass in PHE (kg)
nn_s = 4 # number of coolant nodes in PHE
mn_s = m_s / nn_s # coolant mass per node (kg)
scp_s = 2.39E-3 # cp_s/m_s; % specific heat capacity of coolant (MJ/(kg-C) ORNL-TM-0728 p.8
A_phe = 2.359E+01 # effective area for heat transfer (primary and secondary, m^2) ORNL-TM-0728 p.164
ha_p = 6.480E-01 # heat transfer*area coefficient from primary to tubes (MW/C) ORNL-TM-1647 p.3
ha_s = 3.060E-01 # heat transfer*area coefficient from tubes to secondary (MW/C) ORNL-TM-1647 p.3
# Primary Side
mcp_pn = mn_p * cp_p # (mass of material x heat capacity of material) of fuel salt per lump in MW-s/°C
hA_pn = ha_p / nn_s # 3.030; % (primary to tube heat transfer coeff x heat transfer area) in MW/°C
# Tubes - DONE
nn_t = 2 # number of nodes of tubes in the model
rho_tube = 8.7745E+03 # (kg/m^3) density of INOR-8 ORNL-TM-0728 p.20
m_tn = (v_tube - v_cool) * in_m * rho_tube / nn_t # mass of tubes (kg)
scp_t = 5.778E-04 # specific heat capacity of tubes (MJ/(kg-C)) ORNL-TM-0728 p.20
mcp_tn = m_tn * scp_t # mass*(heat capacity) of tubes per lump in MW-s/°C
# Secondary Side - DONE
mcp_sn = mn_s * scp_s # (mass of material x heat capacity of material) of coolant salt per lump in MW-s/°C
hA_sn = ha_s / nn_s # (tube to secondary heat transfer coeff x heat transfer area) in MW/°C
# Initial conditions - DONE
# Primary nodes
Tp_in = T0_f2 # in °C ORNL-TM-1647 p.2
T0_p4 = Tf_in # 6.5444E+02; % in °C 6.461904761904777e+02; ORNL-TM-1647 p.2
T0_p1 = Tp_in + (T0_p4 - Tp_in) / 4 # in °C
T0_p2 = Tp_in + 2 * (T0_p4 - Tp_in) / 4 # in °C
T0_p3 = Tp_in + 3 * (T0_p4 - Tp_in) / 4 # in °C
# Secondary nodes
Ts_in = 5.4611E+02 # in °C ORNL-TM-1647 p.2
T0_s4 = 5.7939E+02 # in °C ORNL-TM-1647 p.2
T0_s1 = Ts_in + (T0_s4 - Ts_in) / nn_s # in °C
T0_s2 = Ts_in + 2 * (T0_s4 - Ts_in) / nn_s # in °C
T0_s3 = Ts_in + 3 * (T0_s4 - Ts_in) / nn_s # in °C
# Tube nodes
T0_t1 = (T0_p1 * hA_pn + T0_s3 * hA_sn) / (hA_pn + hA_sn) # in °C
T0_t2 = (T0_p3 * hA_pn + T0_s1 * hA_sn) / (hA_pn + hA_sn) # in °C
# Radiator Parameters - DONE
# Initial conditions - DONE
# Primary nodes
Trp_in = T0_s4 # 5.933E+02; % in °C ORNL-TM-1647 p.2
T0_rp = Ts_in # in °C ORNL-TM-1647 p.2
# Secondary nodes - DONE
Trs_in = 37.78 # (C) air inlet temperature ORNL-TM-1647 p.2
T0_rs = 148.9 # (C) air exit temperature ORNL-TM-1647 p.2
# Radiator Geometry
od_rad = 0.01905 # (m) outer diameter of tubes in the radiator ORNL-TM-0728 p.296
tube_wall_thick = 0.0018288 # (m) thickness of tubes in the radiator ORNL-TM-0728 p.296
id_rad = od_rad - 2 * tube_wall_thick
n_rtubes = 120 # number of tubes in the radiator (rows times tubes per row) ORNL-TM-0728 p.296
l_rtube = 9.144 # (m) length of tubes in the radiator ORNL-TM-0728 p.296
v_rp = pi * (id_rad / 2) ** 2 * l_rtube * n_rtubes # volume available to salt in the radiator
# v_rtube = pi * (od_rad / 2) ** 2 * l_rtube * n_rtubes - v_rp # volume of metal in radiator tubes *TUBES NOT MODELED
n_tpr = 12 # number of tubes per row in the radiator matrix
n_row = 10 # number rows in the radiator matrix
tube_space = 0.0381 # (m) spacing between tubes and rows of matrix
v_rs = (n_row * od_rad + (n_row - 1) * tube_space) * (n_tpr * od_rad + (n_tpr - 1) * tube_space) * l_rtube # volume of air inside radiator
# PRIMARY FLOW PARAMETERS - DONE
W_rp = W_s # coolant salt flow rate (kg/s)
m_rp = v_rp * rho_s # coolant salt mass in rad (kg)
nn_rp = 1 # number of coolant salt nodes in the radiator
mn_rp = m_rp / nn_rp # coolant mass per node (kg)
cp_rp = scp_s # coolant specific heat capacity (MJ/(kg-C))
# SECONDARY FLOW PARAMETERS - DONE
vdot_rs = 94.389 # ORNL-TM-0728 p. 296; 78.82; % air volume flow rate (m^3/s) ORNL-TM-1647 p.2
rho_rs = 1.1237 # air density (kg/m^3) REFPROP (310K and 0.1MPa)
W_rs = vdot_rs * rho_rs # air flow rate (kg/s)
m_rs = v_rs * rho_rs # coolant air mass in rad (kg)
nn_rs = 1 # number of coolant nodes in rad
mn_rs = m_rs / nn_rs # coolant mass per node (kg)
scp_rs = 1.0085E-3 # (MJ/kg-C) specific heat capacity of air at (air_out+air_in)/2 REFPROP
A_rad = 6.503E1 # (m^2) surface area of the radiator ORNL-TM-0728 p.14
h_roverall = P / A_rad / ((T0_rp + Trp_in) / 2 - (T0_rs + Trs_in) / 2) # cald as: P/A_rad/((T0_rp+Trp_in)/2-(T0_rs+Trs_in)/2) 3.168E-4; % (MW/m^2-C) polimi thesis
# Primary Side
mcp_rpn = mn_rp * cp_rp # (mass of material x heat capacity of material) of fuel salt per lump in MW-s/°C
hA_rpn = h_roverall * A_rad / nn_rs # 3.030; % (primary to secondary heat transfer coeff x heat transfer area) in MW/°C
# Secondary Side - DONE
mcp_rsn = mn_rs * scp_rs # (mass of material x heat capacity of material) of coolant salt per lump in MW-s/°C
hA_rsn = h_roverall * A_rad / nn_rs # (tube to secondary heat transfer coeff x heat transfer area) in MW/°C
# Pure time delays between components - DONE
tau_hx_c = 8.67 # (sec) delay from hx to core TDAMSRE p.6
tau_c_hx = 3.77 # (sec) subtracted 1 sec for external loop power generation node resident time; delay from core to fuel hx TDAMSRE p.6
tau_hx_r = 4.71 # (sec) fertile hx to core TDAMSRE p.6
tau_r_hx = 8.24 # (sec) core to fertile hx TDAMSRE p.6
first_val = (rho_0 - beta_t) * n_frac0 / Lam + lam[0] * C0[0] + lam[1] * C0[1] + lam[2] * C0[2] + lam[3] * C0[3] + lam[4] * C0[4] + lam[5] * C0[5]

94
dynamic_model/plot_results.py Normal file
View file

@ -0,0 +1,94 @@
import matplotlib.pyplot as plt
from parameters import *
# 1: radiatior coolant inlet temp
# 2: radiator coolant outlet temp
# 3: radiator air outlet temp
# 4: heat exchanger fuel inlet temp
# 5: heat exchanger fuel node 1 temp
# 6: heat exchanger fuel node 2 temp
# 7: heat exchanger fuel node 3 temp
# 8: heat exchanger fuel node 4 temp
# 9: heat exchanger tube node 1 temp
# 10: heat exchanger tube node 2 temp
# 11: heat exchanger coolant inlet temp
# 12: heat exchanger coolant node 1
# 13: heat exchanger coolant node 2
# 14: heat exchanger coolant node 3
# 15: heat exchanger coolant node 4
# 16: core fuel inlet temp
# 17: k
# 18: C1
# 19: C2
# 20: C3
# 21: C4
# 22: C5
# 23: C6
# 24: core graphite temp
# 25: core fuel node 1 temp
# 26: core fuel node 2 temp
filename1 = "sim_out_1000.0_1.txt"
filename5 = "sim_out_1000.0_5.txt"
filename8 = "sim_out_1000.0_8.txt"
sol1 = []
k_file = open(filename1, 'r')
for k in k_file.readlines():
sol1.append(k.split())
k_file.close()
sol5 = []
k_file = open(filename5, 'r')
for k in k_file.readlines():
sol5.append(k.split())
k_file.close()
sol8 = []
k_file = open(filename8, 'r')
for k in k_file.readlines():
sol8.append(k.split())
k_file.close()
sol1 = [[float(j) for j in s] for s in sol1]
sol5 = [[float(j) for j in s] for s in sol5]
sol8 = [[float(j) for j in s] for s in sol8]
k = 10000
test_pow1 = [(1*s[13]-1) for s in sol1]
test_pow5 = [(5*s[13]-5) for s in sol5]
test_pow8 = [(8*s[13]-8) for s in sol8]
tidx1 = [t[0] for t in enumerate(sol1) if (t[1][0] >= 500.00 and t[1][0] <= 800.00)]
tidx5 = [t[0] for t in enumerate(sol5) if (t[1][0] >= 500.00 and t[1][0] <= 800.00)]
tidx8 = [t[0] for t in enumerate(sol8) if (t[1][0] >= 500.00 and t[1][0] <= 800.00)]
t1 = [s[0] for s in sol1[tidx1[0]:tidx1[-1]]]
t5 = [s[0] for s in sol5[tidx5[0]:tidx5[-1]]]
t8 = [s[0] for s in sol8[tidx8[0]:tidx8[-1]]]
# Create a figure and a 3x1 grid of subplots
fig, axs = plt.subplots(3, 1, figsize=(6, 12))
# Plot data on the first subplot
axs[0].plot(t1, test_pow1[tidx1[0]:tidx1[-1]])
axs[0].set_title('1 MW')
#axs[0].set_xlabel('x')
axs[0].set_ylabel('dP')
axs[0].set_xticklabels([])
axs[0].set_xticks([])
# Plot data on the second subplot
axs[1].plot(t5, test_pow5[tidx5[0]:tidx5[-1]])
axs[1].set_title('5 MW')
#axs[1].set_xlabel('x')
axs[1].set_ylabel('dP')
axs[1].set_xticklabels([])
axs[1].set_xticks([])
# Plot data on the third subplot
axs[2].plot(t8, test_pow8[tidx8[0]:tidx8[-1]])
axs[2].set_title('8 MW')
axs[2].set_xlabel('x')
axs[2].set_ylabel('dP')
plt.tight_layout()
plt.show()