API tools faq
paste
casualguitar

Code for Perse95

casualguitar
May 9th, 2022 (edited)
103
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 8.90 KB | None | 0 0
raw download clone embed print report
  1. import numpy as np
  2. from scipy.integrate import solve_ivp
  3. import matplotlib.pyplot as plt
  4. import math
  5.  
  6. def initialize_system(*constants):
  7.     ##sets the appropriate initial conditions
  8.     #which is actually just a vector of constants (initially all nodes have the same initial condition values)
  9.  
  10.     n, dz, y0CO2, y0H2O, Tg0, Tb, MCO2, MH2O = constants #unpack constants
  11.  
  12.     #Initial conditions values are the same at all positions except boundary/inlet at t=0
  13.  
  14.     x0 = np.asarray([[y0CO2, y0H2O, Tg0, Tb, MCO2, MH2O]]*n) #2D array of initial condition values (setting initial conditions at all values of n equal to the same values)
  15.  
  16.     #Initially, all initial values for internal nodes are the same at all points
  17.     #i.e.
  18.     # y0CO2 = y1CO2 = ... = ynCO2
  19.     # y0H2O = y1H20 = ... = ynH20
  20.     # etc
  21.  
  22.     #x0 is a 2D array that looks like this:
  23.     # n=0: |y0CO2, y0H2O, Tg0, Tb0, M0CO2, M1H2O|
  24.     # n=1: |y1CO2, y1H2O, Tg1, Tb1, M1CO2, M2H2O|
  25.     # n=2: |y2CO2, y2H2O, Tg2, Tb2, M2CO2, M3H2O|
  26.     # ...   ..... ......  ...  ...  .....  .....
  27.     # n=n: |ynCO2, ynH2O, Tgn, Tbn, MnCO2, MnH2O|
  28.  
  29.     #It is then reshaped to be a column vector in the form:
  30.     # [y0CO2, y0H2O, Tg0, Tb0, M0CO2, M1H2O,y1CO2, y1H2O, Tg1, Tb1, M1CO2, M2H2O, y2CO2, y2H2O, Tg2, Tb2, M2CO2, M3H2O,...ynCO2, ynH2O, Tgn, Tbn, MnCO2, MnH2O ]
  31.  
  32.     # Meaning that we do not need the for loop setup in initialize_system (possibly wrong?)
  33.     # because we have constants everywhere initially. So we dont need the for loop below unless we have different values everywhere initially?
  34.  
  35.  
  36.     # #Initial condition setup (non boundary points)
  37.     # for j in range(1, n):
  38.     #     # CO2 mole fraction (dyi/dt)
  39.     #     x0[j][0] = (mdot*(y0CO2 - x0[1][j]) - MDR_CO2*A_s + x0[1][j]*1)/m
  40.     #     # H2O mole fraction (dyi/dt)
  41.     #     x0[j][1] = (mdot*(y0H2O - x0[1][j]) - MDR_CO2*A_s + x0[1][j]*1)/m
  42.     #     # CO2 solid deposition (dMi/dt)
  43.     #     x0[j][2] = MDR_CO2*A_s/(A_c*dz)
  44.     #     # H2O solid deposition (or liq) (dMi/dt)
  45.     #     x0[j][3] = MDR_H2O*A_s/(A_c*dz)
  46.     #     # Gas phase heat balance (dTg/dt)
  47.     #     x0[j][4] = (mdot*Cp(Tg0 - x0[4][j]) - Q*A_s)/(m*Cp)
  48.     #     # Bed heat balance (dTb/dt)
  49.     #     x0[j][5] = (Q*A_s)/(M_s*Cp)
  50.     #     #mj
  51.     #     mdot = mdot + (rho/x0[4][j])*A_c*dz*x0[4][j] - (MDR_CO2*A_s)
  52.  
  53.     #reshapes 2D array as a1,b1,c1,d1,e1,a2,b2,c2,d2,e2 or a1,...an,b1...bn, c1...cn,
  54.     return x0.reshape(-1)
  55.  
  56.  
  57. def df(t,x,*constants):
  58.     ##the time derivative of x
  59.     dz, number_of_species, number_of_nodes, A_s, Cp = constants  ## unpack constants
  60.  
  61.     ##setting up solution array
  62.     x = x.reshape(number_of_nodes,6)
  63.     dxdt = np.zeros_like(x)
  64.  
  65.     ##Set up boundary conditions here (n=0, t=t i.e. at position 0 at all times)
  66.  
  67.     ##BCs here
  68.     ##boundary condition values at all times (inlet values)
  69.     mjin = 1 ##molar flow in
  70.     TgIn = 1 ##temperature of flow in
  71.     yInCO2 = 0.1 ##mole fraction of N2
  72.     yInH2O = 0.1 ##mole fraction of H2O
  73.     ##density of components (not needed)
  74.     Cp = 1 ##heat capacity (constant for now)
  75.  
  76.     #boundary ODE setup
  77.     # setting up ODE for inlet position (n=0, t=t). This uses the inlet values, hence
  78.     # the need for both boundary and initial condition setup
  79.  
  80.     # note: the reason I set the boundary conditions here rather than in initialize_system
  81.     # is because I need to reference the solution array to set up the ODEs (initialize_system does not take
  82.     # x as an argument meaning I could not do this, but df does so I set them up there)
  83.  
  84.     #boundary/inlet ODEs
  85.     # CO2 mole fraction (dyi/dt)
  86.     dxdt[0][0] = (mjin * (yInCO2 - x[1][0]) - MDR_CO2 * A_s + x[1][0] * 1) / m
  87.     # H2O mole fraction (dyi/dt)
  88.     dxdt[0][1] = (mjin * (yInH2O - x[1][0]) - MDR_H2O * A_s + x[1][0] * 1) / m
  89.     # CO2 solid deposition (dMi/dt)
  90.     dxdt[0][2] = MDR_CO2 * A_s / (A_c * dz)
  91.     # H2O solid deposition (dMi/dt)
  92.     dxdt[0][3] = MDR_H2O * A_s / (A_c * dz)
  93.     # Gas phase heat balance (dTg/dt)
  94.     dxdt[0][4] = (mjin * Cp(TgIn - x0[4][0]) - Q * A_s) / (m * Cp)
  95.     # Bed heat balance (dTb/dt)
  96.     dxdt[0][5] = (Q * A_s) / (M_s * Cp)
  97.     # mj
  98.     mdot = mjin + (rho / x[4][0]) * A_c * dz * dxdt[0][4] - (MDR_CO2 * A_s)
  99.  
  100.     ## internal nodes (and end node which is the same as the internal nodes)
  101.     for j in range(1,n): ##start at index 1 (to avoid BCs), go to n (last node)
  102.         ## CO2 mole fraction (dyi/dt)
  103.         dxdt[j][0] = (mdot*(x[j-1][0] - x[j][0]) - MDR_CO2*A_s + x[j][0]*1)/m
  104.         ## H2O mole fraction (dyi/dt)
  105.         dxdt[j][1] = (mdot*(x[j-1][1] - x[j][1]) - MDR_CO2*A_s + x[j][1]*1)/m
  106.         ## CO2 solid deposition (dMi/dt)
  107.         dxdt[j][2] = MDR_CO2*A_s/(A_c*dz)
  108.         ## H2O solid deposition (or liq) (dMi/dt)
  109.         dxdt[j][3] = MDR_H2O*A_s/(A_c*dz)
  110.         ## Gas phase heat balance (dTg/dt)
  111.         dxdt[j][4] = (mdot*Cp(x[j-1][4] - x[j][4]) - Q*A_s)/(m*Cp)
  112.         ## Bed heat balance (dTb/dt)
  113.         dxdt[j][5] = (Q*A_s)/(M_s*Cp)
  114.         ## mj
  115.         mdot = mdot + (rho/x[j][4])*A_c*dz*dxdt[j][4] - (MDR_CO2*A_s)
  116.  
  117.     return dxdt.reshape(-1)
  118.  
  119. y0CO2 = 0.1
  120. y0H2O = 0.1
  121. Tg0 = 0.1
  122. Tb = 0.1
  123. MCO2 = 0.1
  124. MH2O = 0.1
  125.  
  126. constants = y0CO2,y0H2O, Tg0, Tb, MCO2, MH2O
  127.  
  128. ##other values here (all this can go as *constants for now but will all be correlations/expressions)
  129. MDR_CO2 = 1 #molar deposition rate
  130. MDR_H2O = 1 #molar deposition rate
  131. m = 1 #will be a correlation for rho and will solve for molar holdup from this
  132. A_c = 1 #cross sectional area
  133. dz = 1 #delta z
  134. Cp = 1 #will be a correlation
  135. Q = 1 #this is a function of T_g, T_b and the heat transfer coefficients U_g and U_b
  136. M_s = 1
  137. rho = 0 #will be a correlation
  138. n = 10
  139. A_s = 1
  140. epsilon = 0.3
  141.  
  142. #create your initial conditions
  143. t0 = 0
  144. tf = 100 #how long you want to integrate for
  145. x0 = initialize_system(constants)
  146.  
  147. #solve the system
  148. sol = solve_ivp(df, [t0,tf], x0, dense_output = True) #dense_output gives you an interpolant (3rd/4th order polynomial)
  149.  
  150. #decide at which times you'd like to check your system
  151. t_eval = [t_i, t_i+1,...t_f]
  152.  
  153. x_eval = [sol.sol(t) for t in t_eval]
  154.  
  155.  
  156.  
  157. #These are the expresstions (currently mostly empty) that will replace some of the constants)
  158. def Ug():
  159.     # U_g = (2.19Re^1/3 + 0.78*Re^0.619)*Sc^1/3 * k_g
  160.     return 0
  161.  
  162. def Ub():
  163.     # U_b = (2.19Re^1/3 + 0.78*Re^0.619)*Sc^1/3 * k_g
  164.     return 0
  165.  
  166. def QgI():
  167.     # QgI= -(U_g*qstar)/((U_g + U_b) + Ustar(T_g-T_b)
  168.     return 0
  169.  
  170. def QIb():
  171.     # QIb = (U_b*qstar)/((U_g + U_b) + Ustar(T_g-T_b)
  172.     return 0
  173.  
  174. def qstar():
  175.     #qstar = M_CO2*deltaH_CO2 + M_H2O*deltaH_H2O
  176.     return 0
  177.  
  178. def Ustar():
  179.     #Ustar = 1/((1/U_g)+ (1/U_b))
  180.     return 0
  181.  
  182. def m():
  183.     #gas phase density * element volume
  184.     return 0
  185.  
  186. def rho_CO2():
  187.     #gas phase density
  188.     #could use thermo here
  189.     return 0
  190.  
  191. def rho_H2O():
  192.     #gas phase density
  193.     #could use thermo here
  194.     return 0
  195.  
  196. def molar_deposition_rate_CO2():
  197.     #
  198.     return 0
  199.  
  200. def molar_deposition_rate_H2O():
  201.     return 0
  202.  
  203. def cp_CO2():
  204.     #thermo
  205.     return 0
  206.  
  207. def cp_H2O():
  208.     #thermo
  209.     return 0
  210.  
  211. def Pr():
  212.     # Pr = mu_g*Cp/k_g
  213.     return 0
  214.  
  215. def Re(m_molar,T, A_c, epsilon,d):
  216.     #Re = 6G0/(a*mu_g*Psi)
  217.     Psi = 1
  218.     g0 = m_molar/A_c
  219.     sa_specific = (1-epsilon)*6/d #surface area per unit volume
  220.     mu = 1 #correlation here #Pascal seconds #Pa.s
  221.     re = 6*g0/(sa_specific*mu)
  222.  
  223.     return re
  224.  
  225. def Nu():
  226.     #Nu = h_fs*d_p/k_g
  227.     return 0
  228.  
  229. def Sh():
  230.     #Sh = k_g*D_p/((1-e_g)Di*Psi)
  231.     #D_i is diffusion coefficient of component i
  232.     return 0
  233.  
  234. def Sc():
  235.     #Sc = mu/rho*D_i
  236.  
  237.  
  238.  
  239.     return 0
  240.  
  241. def k():
  242.     #k_f = (2.19Re^1/3 + 0.78*Re^0.619)*Sc^1/3 * k_g
  243.     return 0
  244.  
  245. def psubCO2(T):
  246.     # sublimation pressure equation
  247.     # only applies up to the triple point pressure Pt
  248.  
  249.     # arguments: Temperature
  250.     Tt = 216.592  # triple point temperature (K)
  251.     Pt = 0.51795  # triple point pressure (MPa)
  252.     a1 = -14.740846  # constant a1
  253.     a2 = 2.4327015  # constant a2
  254.     a3 = -5.3061778  # constant a3
  255.  
  256.     if T>Tt:
  257.         return 0
  258.     else:
  259.         reduced_temperature = T / Tt
  260.         term1 = a1 * (1 - reduced_temperature)
  261.         term2 = a2 * (1 - reduced_temperature) ** 1.9
  262.         term3 = a3 * (1 - reduced_temperature) ** 2.9
  263.         term4 = (1 / reduced_temperature) * (term1 + term2 + term3)
  264.         term5 = (math.e) ** term4
  265.         psub = term5 * Pt
  266.  
  267.     return psub
  268.  
  269. def pliqH2O(T):
  270.     #liquefaction pressure equation
  271.     #arguments: Temperature
  272.     Pt = 611.657 #Pa
  273.     Tt = 273.16 #K
  274.     a1 = 0.119539337*10**7 #constant a1
  275.     a2 = 0.808183159*10**5 #constant a2
  276.     a3 = 0.333826860*10**4 #constant a3
  277.     b1 = 0.300*10**1 #constant b1
  278.     b2 = 0.2575*10**2 #constant b2
  279.     b3 = 0.103750*10**3 #constant b3
  280.     theta = T/Tt #reduced temperature
  281.  
  282.     term1 = a1*(1-theta**b1)
  283.     term2 = a2*(1-theta**b2)
  284.     term3 = a3*(1-theta**b3)
  285.     term4 = 1 + term1 + term2 + term3
  286.     pLiq = term4*Pt
  287.  
  288.     return pLiq
  289.  
  290.  
  291.  
  292.  
  293.  
  294.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!