API tools faq
paste
casualguitar

numblsoda script

casualguitar
Feb 28th, 2022
175
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 12.48 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. from chemicals import normalize
  5. from thermo import ChemicalConstantsPackage, CEOSGas, CEOSLiquid, PRMIXTranslatedConsistent, FlashVL
  6. from thermo.interaction_parameters import IPDB
  7. import math
  8. from thermo.chemical import Mixture
  9.  
  10. def temperature(H,P):
  11.     # returns temperature given enthalpy (after processing function)
  12.     T = flasher.flash(H=H, P=P, zs=zs, retry=True).T #K
  13.     return T
  14.  
  15. def mass(H, P):
  16.     # returns mass holdup in mol (mass in an element n)
  17.     m = flasher.flash(H=H, P=P, zs=zs, retry=True).rho()*V #mol
  18.     return m
  19.  
  20. def dpdh(H, P):
  21.     #returns the derivative of density w.r.t enthalpy at constant pressure
  22.     res = flasher.flash(H=H, P=P, zs=zs, retry=True)
  23.     if res.phase_count == 1:
  24.         if res.phase == 'L':
  25.             drho_dTf = res.liquid0.drho_dT()
  26.         else:
  27.             drho_dTf = res.gas.drho_dT()
  28.     else:
  29.         drho_dTf = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='P')
  30.     dpdh = drho_dTf/res.dH_dT_P()
  31.     return dpdh
  32.  
  33. def U(H,P,m):
  34.     # Given T, P, m_dot returns U W/m2.K
  35.     air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
  36.     mu = air.mu*1000/MW_AIR #mol/m.s
  37.     cp = air.Cpm #J/mol.K
  38.     kg = air.k #W/m.K
  39.     g0 = m/AREA_BED #mol/m2.s
  40.     a = sa*n/V_TOTAL #m^2/m^3 #QUESTIONABLE
  41.     psi = 1
  42.     beta = 10
  43.     pr = (mu*cp)/kg
  44.     re = (6*g0)/(a*mu*psi)
  45.     hfs = ((2.19*(re**1/3)) + (0.78*(re**0.619)))*(pr**1/3)*(kg)/DIAMETER_PARTICLE
  46.  
  47.     h = 1/((1/hfs) + ((DIAMETER_PARTICLE/beta)/K_ALUMINA))
  48.  
  49.     return h
  50.  
  51. ##Setup of the fluid mixture properties (air)
  52. constants, properties = ChemicalConstantsPackage.from_IDs(['nitrogen', 'oxygen', 'argon'])
  53. kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
  54. eos_kwargs = {'Pcs': constants.Pcs, 'Tcs': constants.Tcs, 'omegas': constants.omegas, 'kijs': kijs}
  55. gas = CEOSGas(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
  56. liquid = CEOSLiquid(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
  57. flasher = FlashVL(constants, properties, liquid=liquid, gas=gas)
  58. zs = normalize([78.08, 20.95, .93])
  59. ################################################################
  60.  
  61. #BED
  62. NUMBER_BEDS = 10 #Physical packed beds
  63. DIAMETER_BED = 0.0779 #metres
  64. HEIGHT_BED = 0.27 #metres
  65. AREA_BED = ((math.pi)/4)*DIAMETER_BED**2 #m^2
  66. VOLUME_BED = AREA_BED*HEIGHT_BED #m^3
  67. V_TOTAL = VOLUME_BED * NUMBER_BEDS #total beds volume (m^3)
  68.  
  69. #PARTICLE
  70. DIAMETER_PARTICLE = 0.005 #metres
  71. VOLUME_PARTICLE = (1/6)*math.pi*(DIAMETER_PARTICLE**3) #particle volume m^3
  72. PARTICLE_NUMBER = math.floor(V_TOTAL/VOLUME_PARTICLE) #total number of particles
  73. SA_PARTICLE = 4*math.pi*((DIAMETER_PARTICLE/2)**2) #m^2
  74. SA_TOTAL = SA_PARTICLE*PARTICLE_NUMBER #total surface area (m^2)
  75. RHO_ALUMINA = 3950 #kg/m3
  76. M_ALUMINA_KG = VOLUME_PARTICLE*PARTICLE_NUMBER*RHO_ALUMINA # total mass of alumina (kg)
  77. MW_ALUMINA = 101.96 #g/mol
  78. M_ALUMINA_MOL = M_ALUMINA_KG*1000/MW_ALUMINA #total mass of alumina (mol)
  79. CP_ALUMINA = 89.72 #J/mol.K
  80. K_ALUMINA = 5 #W/m.K
  81.  
  82. #FLUID
  83. MW_AIR = 29 #g/mol
  84.  
  85. #INITIAL CONDITIONS
  86. T0 = 80  # initial fluid temperature (K)
  87. Ts0 = 80  # initial solid temperature (K)
  88. P = 3000000  # Pa
  89. H0 = flasher.flash(T=T0, P=P, zs=zs, retry = True).H() # bed fluid starting enthalpy (J/mol)
  90.  
  91.  
  92. #BOUNDARY CONDITIONS
  93. Tin = 300  # Inlet fluid temperature (K)
  94. minkg = 0.022 #kg/s
  95. min = minkg*1000/MW_AIR #mol/s
  96. hIn = flasher.flash(T=Tin, P=P, zs=zs, retry = True).H()  # Inlet fluid enthalpy (J/mol)
  97.  
  98. #SPATIAL DEPENDENCE OF QUANTITIES
  99. n = 3 # number of elements in spatial array
  100. V = V_TOTAL/n
  101. sa = SA_TOTAL/n #surface area per unit volume
  102. mAlumina = M_ALUMINA_MOL/n #alumina mass per unit volume
  103.  
  104. # time steps
  105. TIME = 1000
  106. t = (0, TIME)  # time range
  107. # t_eval = np.linspace(0, TIME, 1000 * TIME)  # optional
  108.  
  109.  
  110. def my_system(t, y, n, hIn, min, mAlumina, CP_ALUMINA, sa, V):
  111.     dydt = np.zeros(2 * n) #setting up zeros array
  112.  
  113.     i=0
  114.  
  115.     ## Using thermo
  116.     T = temperature(y[i],P) #initial fluid temperature
  117.     m = mass(y[i],P) #initial mass holdup
  118.  
  119.     #initial boundary values
  120.     dydt[i] = (min * (hIn - y[i]) + (U(y[i],P,min) * sa * (y[i + 1] - T))) / m # dH/dt (eq. 2)
  121.     dydt[i + 1] = -(U(y[i],P,min) * sa * (y[i + 1] - T)) / (mAlumina * CP_ALUMINA) # dTs/dt from eq.3
  122.     dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt eq. 4b
  123.     mj = min - dmdt # mass flow out (eq.4a)
  124.  
  125.     #Solution array y will look like [h0,Ts0,mdot0, h1, Ts1, mdot1,...hn,Tsn,mdotn] at all values of t
  126.     for i in range(2, 2 * n, 2): #starting at index 3, and incrementing by 3 because we are solving for 'triplets' [h,Ts,m] in each loop
  127.  
  128.         ## Using thermo
  129.         T = temperature(y[i],P)
  130.         m = mass(y[i],P)
  131.  
  132.         # [h, TS, mdot]
  133.         dydt[i] = (mj * (y[i - 2] - y[i]) + (U(y[i], P, mj) * sa * (y[i + 1] - T))) / m # dH/dt (eq.2)
  134.         dydt[i + 1] = -(U(y[i], P, mj) * sa * (y[i + 1] - T)) / (mAlumina * CP_ALUMINA) # dTs/dt eq. (3)
  135.         dmdt = dydt[i] * dpdh(y[i], P) * V # Equation 4b
  136.         mj = mj - dmdt # Equation 4a
  137.  
  138.     return dydt
  139.  
  140.  
  141. y0 = np.copy(np.asarray([H0, Ts0] * n)) #setting up initial condition array (H0,Ts0 and m at all points
  142. y0[0] = H0 # initial fluid bed enthalpy
  143. y0[1] = Ts0 # initial solid bed temperature
  144.  
  145.  
  146. parameters = (n, hIn, min, mAlumina, CP_ALUMINA, sa, V)
  147.  
  148. solution = solve_ivp(my_system, t, y0, args=parameters, method='Radau')
  149.  
  150. h_arr = solution.y[0::2]
  151. Ts_arr = solution.y[1::2]
  152. import numpy as np
  153. from scipy.integrate import solve_ivp
  154. import matplotlib.pyplot as plt
  155. from chemicals import normalize
  156. from thermo import ChemicalConstantsPackage, CEOSGas, CEOSLiquid, PRMIXTranslatedConsistent, FlashVL
  157. from thermo.interaction_parameters import IPDB
  158. import math
  159. from thermo.chemical import Mixture
  160.  
  161. def temperature(H,P):
  162.     # returns temperature given enthalpy (after processing function)
  163.     T = flasher.flash(H=H, P=P, zs=zs, retry=True).T #K
  164.     return T
  165.  
  166. def mass(H, P):
  167.     # returns mass holdup in mol (mass in an element n)
  168.     m = flasher.flash(H=H, P=P, zs=zs, retry=True).rho()*V #mol
  169.     return m
  170.  
  171. def dpdh(H, P):
  172.     #returns the derivative of density w.r.t enthalpy at constant pressure
  173.     res = flasher.flash(H=H, P=P, zs=zs, retry=True)
  174.     if res.phase_count == 1:
  175.         if res.phase == 'L':
  176.             drho_dTf = res.liquid0.drho_dT()
  177.         else:
  178.             drho_dTf = res.gas.drho_dT()
  179.     else:
  180.         drho_dTf = res.bulk._equilibrium_derivative(of='rho', wrt='T', const='P')
  181.     dpdh = drho_dTf/res.dH_dT_P()
  182.     return dpdh
  183.  
  184. def U(H,P,m):
  185.     # Given T, P, m_dot returns U W/m2.K
  186.     air = Mixture(['nitrogen', 'oxygen'], Vfgs=[0.79, 0.21], H=H, P=P)
  187.     mu = air.mu*1000/MW_AIR #mol/m.s
  188.     cp = air.Cpm #J/mol.K
  189.     kg = air.k #W/m.K
  190.     g0 = m/AREA_BED #mol/m2.s
  191.     a = sa*n/V_TOTAL #m^2/m^3 #QUESTIONABLE
  192.     psi = 1
  193.     beta = 10
  194.     pr = (mu*cp)/kg
  195.     re = (6*g0)/(a*mu*psi)
  196.     hfs = ((2.19*(re**1/3)) + (0.78*(re**0.619)))*(pr**1/3)*(kg)/DIAMETER_PARTICLE
  197.  
  198.     h = 1/((1/hfs) + ((DIAMETER_PARTICLE/beta)/K_ALUMINA))
  199.  
  200.     return h
  201.  
  202. def vapourfraction(H, P):
  203.     # returns mass holdup in mol (mass in an element n)
  204.     x = flasher.flash(H=H, P=P, zs=zs, retry=True).VF #mol
  205.     return x
  206.  
  207. ##Setup of the fluid mixture properties (air)
  208. constants, properties = ChemicalConstantsPackage.from_IDs(['nitrogen', 'oxygen', 'argon'])
  209. kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
  210. eos_kwargs = {'Pcs': constants.Pcs, 'Tcs': constants.Tcs, 'omegas': constants.omegas, 'kijs': kijs}
  211. gas = CEOSGas(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
  212. liquid = CEOSLiquid(PRMIXTranslatedConsistent, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
  213. flasher = FlashVL(constants, properties, liquid=liquid, gas=gas)
  214. zs = normalize([78.08, 20.95, .93])
  215. ################################################################
  216.  
  217. #BED
  218. NUMBER_BEDS = 10 #Physical packed beds
  219. DIAMETER_BED = 0.0779 #metres
  220. HEIGHT_BED = 0.27 #metres
  221. AREA_BED = ((math.pi)/4)*DIAMETER_BED**2 #m^2
  222. VOLUME_BED = AREA_BED*HEIGHT_BED #m^3
  223. V_TOTAL = VOLUME_BED * NUMBER_BEDS #total beds volume (m^3)
  224.  
  225. #PARTICLE
  226. DIAMETER_PARTICLE = 0.005 #metres
  227. VOLUME_PARTICLE = (1/6)*math.pi*(DIAMETER_PARTICLE**3) #particle volume m^3
  228. PARTICLE_NUMBER = math.floor(V_TOTAL/VOLUME_PARTICLE) #total number of particles
  229. SA_PARTICLE = 4*math.pi*((DIAMETER_PARTICLE/2)**2) #m^2
  230. SA_TOTAL = SA_PARTICLE*PARTICLE_NUMBER #total surface area (m^2)
  231. RHO_ALUMINA = 3950 #kg/m3
  232. M_ALUMINA_KG = VOLUME_PARTICLE*PARTICLE_NUMBER*RHO_ALUMINA # total mass of alumina (kg)
  233. MW_ALUMINA = 101.96 #g/mol
  234. M_ALUMINA_MOL = M_ALUMINA_KG*1000/MW_ALUMINA #total mass of alumina (mol)
  235. CP_ALUMINA = 89.72 #J/mol.K
  236. K_ALUMINA = 5 #W/m.K
  237.  
  238. #FLUID
  239. MW_AIR = 29 #g/mol
  240.  
  241. #INITIAL CONDITIONS
  242. T0 = 80  # initial fluid temperature (K)
  243. Ts0 = 80  # initial solid temperature (K)
  244. P = 3000000  # Pa
  245. H0 = flasher.flash(T=T0, P=P, zs=zs, retry = True).H() # bed fluid starting enthalpy (J/mol)
  246.  
  247.  
  248. #BOUNDARY CONDITIONS
  249. Tin = 300  # Inlet fluid temperature (K)
  250. minkg = 0.022 #kg/s
  251. min = minkg*1000/MW_AIR #mol/s
  252. hIn = flasher.flash(T=Tin, P=P, zs=zs, retry = True).H()  # Inlet fluid enthalpy (J/mol)
  253.  
  254. #SPATIAL DEPENDENCE OF QUANTITIES
  255. n = 3 # number of elements in spatial array
  256. V = V_TOTAL/n
  257. sa = SA_TOTAL/n #surface area per unit volume
  258. mAlumina = M_ALUMINA_MOL/n #alumina mass per unit volume
  259.  
  260. # time steps
  261. TIME = 1000
  262. t = (0, TIME)  # time range
  263. # t_eval = np.linspace(0, TIME, 1000 * TIME)  # optional
  264.  
  265.  
  266. def my_system(t, y, n, hIn, min, mAlumina, CP_ALUMINA, sa, V):
  267.     dydt = np.zeros(2 * n) #setting up zeros array
  268.  
  269.     i=0
  270.  
  271.     ## Using thermo
  272.     T = temperature(y[i],P) #initial fluid temperature
  273.     m = mass(y[i],P) #initial mass holdup
  274.  
  275.     #initial boundary values
  276.     dydt[i] = (min * (hIn - y[i]) + (U(y[i],P,min) * sa * (y[i + 1] - T))) / m # dH/dt (eq. 2)
  277.     dydt[i + 1] = -(U(y[i],P,min) * sa * (y[i + 1] - T)) / (mAlumina * CP_ALUMINA) # dTs/dt from eq.3
  278.     dmdt = dydt[i] * dpdh(y[i], P) * V # dm/dt eq. 4b
  279.     mj = min - dmdt # mass flow out (eq.4a)
  280.  
  281.     #Solution array y will look like [h0,Ts0,mdot0, h1, Ts1, mdot1,...hn,Tsn,mdotn] at all values of t
  282.     for i in range(2, 2 * n, 2): #starting at index 3, and incrementing by 3 because we are solving for 'triplets' [h,Ts,m] in each loop
  283.  
  284.         ## Using thermo
  285.         T = temperature(y[i],P)
  286.         m = mass(y[i],P)
  287.  
  288.         # [h, TS, mdot]
  289.         dydt[i] = (mj * (y[i - 2] - y[i]) + (U(y[i], P, mj) * sa * (y[i + 1] - T))) / m # dH/dt (eq.2)
  290.         dydt[i + 1] = -(U(y[i], P, mj) * sa * (y[i + 1] - T)) / (mAlumina * CP_ALUMINA) # dTs/dt eq. (3)
  291.         dmdt = dydt[i] * dpdh(y[i], P) * V # Equation 4b
  292.         mj = mj - dmdt # Equation 4a
  293.  
  294.     return dydt
  295.  
  296.  
  297. y0 = np.copy(np.asarray([H0, Ts0] * n)) #setting up initial condition array (H0,Ts0 and m at all points
  298. y0[0] = H0 # initial fluid bed enthalpy
  299. y0[1] = Ts0 # initial solid bed temperature
  300.  
  301.  
  302. parameters = (n, hIn, min, mAlumina, CP_ALUMINA, sa, V)
  303.  
  304. solution = solve_ivp(my_system, t, y0, args=parameters, method='LSODA')
  305.  
  306. #AFTERPROCESSING
  307. h_arr = solution.y[0::2]
  308. Ts_arr = solution.y[1::2]
  309.  
  310. time = solution.t
  311. # dt = 2*(time[2:]-time[:-2])
  312.  
  313. dydt = np.zeros((2 * n, time.size))
  314. for i in range(time.size):
  315.     dydt[:, i] = my_system(t, solution.y[:, i], n, hIn, min, mAlumina, CP_ALUMINA, sa, V)
  316.  
  317. dhdt_arr = dydt[0::2]
  318. dTsdt_arr = dydt[1::2]
  319.  
  320. dmdt_arr = np.zeros_like(dhdt_arr) #change in mass holdup
  321. T_arr = np.zeros_like(dhdt_arr) # fluid temperature
  322. m_hold_arr = np.zeros_like(dhdt_arr) # mass hold up
  323. U_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
  324. mj_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
  325. x_arr = np.zeros_like(dhdt_arr)
  326.  
  327. for i in range(time.shape[0]): #iteratung through 2D array time x #positions
  328.     for j in range(n): #yep
  329.  
  330.         T = temperature(h_arr[j,i],P)
  331.  
  332.         dmdt_arr[j, i] = dhdt_arr[j, i] * dpdh(h_arr[j, i], P) * V
  333.         T_arr[j,i] = T
  334.         m_hold_arr[j, i] = mass(h_arr[j, i], P)
  335.         U_arr[j,i] = mAlumina*CP_ALUMINA*dTsdt_arr[j,i]/(sa*(T_arr[j,i] - Ts_arr[j,i]))
  336.         #x_arr[j,i] = vapourfraction(h_arr[j,i], P)
  337.  
  338.  
  339. # creating Tf_arr from Ts_arr
  340. Tf_arr = np.zeros_like(dhdt_arr)
  341. Tf_arr = Ts_arr + dTsdt_arr * (mAlumina * CP_ALUMINA) / (U_arr * sa)
  342. # Tf_arr[:, 0] = Ts0  # I.C cannot be calculated as there is no derivative at the boundary. Manually entering here
  343.  
  344. print("simulation complete")
  345.  
Advertisement
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!