From 5930aeceec504adfcfcad4237d402ae671865734 Mon Sep 17 00:00:00 2001 From: lukelabrie Date: Mon, 7 Aug 2023 16:52:55 +0200 Subject: [PATCH] add dynamic model --- dynamic_model/main.py | 249 ++++++++++++++++++++++++++++++++++ dynamic_model/parameters.py | 234 ++++++++++++++++++++++++++++++++ dynamic_model/plot_results.py | 94 +++++++++++++ 3 files changed, 577 insertions(+) create mode 100644 dynamic_model/main.py create mode 100644 dynamic_model/parameters.py create mode 100644 dynamic_model/plot_results.py diff --git a/dynamic_model/main.py b/dynamic_model/main.py new file mode 100644 index 0000000..b6f23b1 --- /dev/null +++ b/dynamic_model/main.py @@ -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=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() \ No newline at end of file diff --git a/dynamic_model/parameters.py b/dynamic_model/parameters.py new file mode 100644 index 0000000..365490a --- /dev/null +++ b/dynamic_model/parameters.py @@ -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] \ No newline at end of file diff --git a/dynamic_model/plot_results.py b/dynamic_model/plot_results.py new file mode 100644 index 0000000..b2133dc --- /dev/null +++ b/dynamic_model/plot_results.py @@ -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() \ No newline at end of file