API tools faq
paste
casualguitar

What is 'y' here? Not defined anywhere in the script

casualguitar
Feb 27th, 2022
211
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 12.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. 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.  
  153. time = solution.t
  154. # dt = 2*(time[2:]-time[:-2])
  155.  
  156. dydt = np.zeros((2 * n, time.size))
  157. for i in range(time.size):
  158.     dydt[:, i] = my_system(t, solution.y[:, i], n, hIn, min, mAlumina, CP_ALUMINA, sa, V)
  159.  
  160. dhdt_arr = dydt[0::2]
  161. dTsdt_arr = dydt[1::2]
  162.  
  163. dmdt_arr = np.zeros_like(dhdt_arr) #change in mass holdup
  164. T_arr = np.zeros_like(dhdt_arr) # fluid temperature
  165. m_hold_arr = np.zeros_like(dhdt_arr) # mass hold up
  166. U_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
  167. mj_arr = np.zeros_like(dhdt_arr) # heat transfer coefficient array
  168.  
  169. for i in range(time.shape[0]): #iteratung through 2D array time x #positions
  170.     for j in range(n): #yep
  171.  
  172.         T = temperature(h_arr[j,i],P)
  173.  
  174.         dmdt_arr[j, i] = dhdt_arr[j, i] * dpdh(h_arr[j, i], P) * V
  175.         T_arr[j,i] = T
  176.         m_hold_arr[j, i] = mass(h_arr[j, i], P)
  177.         U_arr[j,i] = mAlumina*CP_ALUMINA*dTsdt_arr[j,i]/(sa*(T_arr[j,i] - Ts_arr[j,i]))
  178.  
  179. # creating Tf_arr from Ts_arr
  180. Tf_arr = np.zeros_like(dhdt_arr)
  181. Tf_arr = Ts_arr + dTsdt_arr * (mAlumina * CP_ALUMINA) / (U_arr * sa)
  182. # Tf_arr[:, 0] = Ts0  # I.C cannot be calculated as there is no derivative at the boundary. Manually entering here
  183.  
  184. print("simulation complete")
  185.  
  186. # plt.figure()
  187. # snapshot = 100  # snapshot time spacing: 1000=1s
  188. # style = ['-', '--', '-.']
  189. # color = ['b','r','g']
  190. # for i in range(len(style)):
  191. #     plt.plot(Tf_arr[:, i*snapshot], style[i],
  192. #              color=color[i],
  193. #              label="t="+str(i*snapshot/1000)+"s")
  194. #
  195. # plt.ylabel(r'Fluid Temperature (K)')
  196. # plt.xlabel(r'position(j)')
  197. # plt.legend(loc='lower left')
  198.  
  199.  
  200. # plt.plot(time,h_arr[0, :], '-.', color='b', label="position 0 h")
  201. # plt.plot(time,h_arr[1, :], '--', color='b', label="position 1")
  202. # plt.plot(time,h_arr[2, :], ':', color='b', label="position 2")
  203. # plt.plot(time,Tf_arr[0, :], '-.', color='g', label="position 0 T")
  204. #plt.plot(time,x_arr[0, :], '-.', color='r', label="tank 1 quality")
  205. #plt.plot(time,x_arr[1, :], '-.', color='g', label="tank 2 quality")
  206. #plt.plot(time,x_arr[2, :], '-.', color='b', label="tank 3 quality")
  207. '''
  208. plt.figure()
  209. plt.plot(time,dmdt_arr[0, :], '-', color='r', label="tank 1 dm/dt")
  210. plt.plot(time,dmdt_arr[1, :], '--', color='g', label="tank 2 dm/dt")
  211. plt.plot(time,dmdt_arr[2, :], '-.', color='b', label="tank 3 dm/dt")
  212. plt.plot(time,dmdt_arr[3, :], '--', color='m', label="tank 4 dm/dt")
  213. plt.plot(time,dmdt_arr[4, :], '-.', color='c', label="tank 5 dm/dt")
  214.  
  215. # plt.plot(time,Tf_arr[1, :], '--', color='g', label="position 1")
  216. # plt.plot(time,Tf_arr[2, :], ':', color='g', label="position 2")
  217. # from plot_together import plot_together
  218. # plot_together(data1[:,[0,1]:], '--', color='b', label="position 0")
  219. # plot_together(data1[:,[0,2],:], '--', color='b', label="position 1")
  220. # plot_together(data1[:,[0,3],:], ':', color='b', label="position 2")
  221. plt.ylabel(r'dm/dt')
  222. plt.xlabel(r'Time')
  223. # plt.ylim((80, 91))
  224. plt.grid('both')
  225. plt.title("dm/dt for all tanks")
  226. plt.legend(loc='lower right')
  227. '''
  228.  
  229. # plt.figure()
  230. # plt.plot(time,mdot_arr[0, :], '-', color='b', label="tank 1 mdot")
  231. # plt.plot(time,mdot_arr[1, :], '--', color='b', label="tank 2 mdot")
  232. # plt.plot(time,mdot_arr[2, :], '-.', color='b', label="tank 3 mdot")
  233. #plt.plot(time,mdot_arr[3, :], '--', color='b', label="tank 4 mdot")
  234. #plt.plot(time,mdot_arr[4, :], '-.', color='b', label="tank 5 mdot")
  235.  
  236. # plt.ylabel(r'Mass Flow (mol/s)')
  237. # plt.xlabel(r'Time')
  238. # plt.ylim((80, 91))
  239. # plt.grid('both')
  240. # plt.title("m dot for all tanks")
  241. # plt.legend(loc='lower right')
  242. #
  243. #########################################################
  244. # plt.figure()
  245. # plt.plot(time,m_arr[0, :], '-', color='b', label="tank 1  m_flow")
  246. # plt.plot(time,m_arr[1, :], '--', color='b', label="tank 2 m_flow")
  247. # plt.plot(time,m_arr[2, :], '-.', color='b', label="tank 3 m_flow")
  248. #plt.plot(time,m_arr[3, :], '--', color='b', label="tank 4 m_flow")
  249. #plt.plot(time,m_arr[4, :], '-.', color='b', label="tank 5 m_flow")
  250.  
  251. # plt.ylabel(r'Mass flow')
  252. # plt.xlabel(r'Time')
  253. # plt.ylim((80, 91))
  254. # plt.grid('both')
  255. # plt.title("Mass flow for all tanks")
  256. # plt.legend(loc='lower right')
  257.  
  258. ###########################################################
  259. plt.figure()
  260. plt.plot(time,T_arr[0, :], '-', color='b' , label="tank 1 Fluid temperature")
  261. plt.plot(time,T_arr[1, :], '--', color='b', label="tank 2 Fluid temperature")
  262. plt.plot(time,T_arr[2, :], '-.', color='b', label="tank 3 Fluid temperature")
  263. plt.plot(time,Tf_arr[0, :], '-', color='r' , label="tank 1 Fluid temperature")
  264. plt.plot(time,Tf_arr[1, :], '--', color='r', label="tank 2 Fluid temperature")
  265. plt.plot(time,Tf_arr[2, :], '-.', color='r', label="tank 3 Fluid temperature")
  266. #plt.plot(time,T_arr[3, :], '--', color='b', label="tank 4 Fluid temperature")
  267. #plt.plot(time,T_arr[4, :], '-.', color='b', label="tank 5 Fluid temperature")
  268.  
  269. plt.ylabel(r'Fluid temperature')
  270. plt.xlabel(r'Time')
  271. # plt.ylim((80, 91))
  272. plt.grid('both')
  273. plt.title("Fluid temperature for all tanks")
  274. plt.legend(loc='lower right')
  275.  
  276. ###########################################################
  277. plt.figure()
  278. plt.plot(time,m_hold_arr[0, :], '-', color='b' , label="tank 1 Mass hold up")
  279. plt.plot(time,m_hold_arr[1, :], '--', color='b', label="tank 2 Mass hold up")
  280. plt.plot(time,m_hold_arr[2, :], '-.', color='b', label="tank 3 Mass hold up")
  281. #plt.plot(time,m_hold_arr[3, :], '--', color='b', label="tank 4 Mass hold up")
  282. #plt.plot(time,m_hold_arr[4, :], '-.', color='b', label="tank 5 Mass hold up")
  283.  
  284. plt.ylabel(r'Mass hold up (mol)')
  285. plt.xlabel(r'Time')
  286. # plt.ylim((80, 91))
  287. plt.grid('both')
  288. plt.title("Mass hold up for all tanks")
  289. plt.legend(loc='lower right')
  290.  
  291. ###########################################################
  292. # plt.figure()
  293. # plt.plot(time,mdot_arr[0, :], '-', color='b' , label="tank 1 Mass flow rate")
  294. # plt.plot(time,mdot_arr[1, :], '--', color='b', label="tank 2 Mass flow rate")
  295. # plt.plot(time,mdot_arr[2, :], '-.', color='b', label="tank 3 Mass flow rate")
  296. #plt.plot(time,mdot_arr[3, :], '--', color='b', label="tank 4 Mass flow rate")
  297. #plt.plot(time,mdot_arr[4, :], '-.', color='b', label="tank 5 Mass flow rate")
  298.  
  299. # plt.ylabel(r'Mass flow rate')
  300. # plt.xlabel(r'Time')
  301. # plt.ylim((80, 91))
  302. # plt.grid('both')
  303. # plt.title("Mass flow rate for all tanks")
  304. # plt.legend(loc='lower right')
  305.  
  306. ###########################################################
  307. plt.figure()
  308. plt.plot(time,Ts_arr[0, :], '-', color='b' , label="tank 1 Solid Temperature")
  309. plt.plot(time,Ts_arr[1, :], '--', color='b', label="tank 2 Solid Temperature")
  310. plt.plot(time,Ts_arr[2, :], '-.', color='b', label="tank 3 Solid Temperature")
  311. #plt.plot(time,Ts_arr[3, :], '--', color='b', label="tank 4 Solid Temperature")
  312. #plt.plot(time,Ts_arr[4, :], '-.', color='b', label="tank 5 Solid Temperature")
  313.  
  314. plt.ylabel(r'Solid Temperature')
  315. plt.xlabel(r'Time')
  316. # plt.ylim((80, 91))
  317. plt.grid('both')
  318. plt.title("Solid Temperature for all tanks")
  319. plt.legend(loc='lower right')
  320.  
  321. ###########################################################
  322. plt.figure()
  323. plt.plot(time,h_arr[0, :], '-', color='b', label="tank 1 h")
  324. plt.plot(time,h_arr[1, :], '--', color='b', label="tank 2 h")
  325. plt.plot(time,h_arr[2, :], '-.', color='b', label="tank 3 h")
  326. #plt.plot(time,h_arr[3, :], '--', color='b', label="tank 4 mdot")
  327. #plt.plot(time,h_arr[4, :], '-.', color='b', label="tank 5 mdot")
  328.  
  329. plt.ylabel(r'h (J/mol)')
  330. plt.xlabel(r'Time')
  331. # plt.ylim((80, 91))
  332. plt.grid('both')
  333. plt.title("h for all tanks")
  334. plt.legend(loc='lower right')
  335.  
  336. plt.show()
  337.  
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!