PYTHON   34
update plot2
Guest on 1st February 2023 01:29:17 AM


  1. #%%
  2. #Libraries
  3. import numpy as np
  4. from scipy.integrate import odeint
  5. import matplotlib.pyplot as plt
  6. import matplotlib as mpl
  7. mpl.rcParams.update({'font.size': 13, 'lines.linewidth': 2.5})
  8. from matplotlib.widgets import Slider, Button
  9.  
  10. #%%
  11. #Explicit equations
  12. Cao=3.99
  13. Cco= 1.94
  14. Cto= 2.09
  15. Vo= 2.38
  16. vo= 6
  17. Cai= 2.38
  18. Nci= 7.76
  19. Nti=7.09
  20. UA=5
  21. def ODEfun(Yfuncvec,t,Cao,Cco,Cto,Vo,vo,Cai,Nci,Nti,UA):
  22.     Ca= Yfuncvec[0]
  23.     Nc= Yfuncvec[1]
  24.     Nt= Yfuncvec[2]
  25.     Can= Yfuncvec[3]
  26.     T= Yfuncvec[4]
  27.     Cas= Yfuncvec[5]
  28.     Ncs= Yfuncvec[6]
  29.     Nts= Yfuncvec[7]
  30.     Cans= Yfuncvec[8]
  31.     Ts= Yfuncvec[9]
  32.       #Explicit Equation Inline
  33.     V = Vo + (vo*10**-5*t);
  34.     Na= Ca*V;
  35.     Nan= Can*V;
  36.     Fao = Cao*vo*10**-5;
  37.     Fco = Cco*vo*10**-5;
  38.     Fto = Cto*vo*10**-5;
  39.     Cpa = 245.75;
  40.     Cpc = 161.3;
  41.     Cpt = 165.60;
  42.     Cpan = 231;
  43.     k= (4.01*10**3)*np.exp(-29000/(8.314*T))
  44.     ra= -k*(Ca)**2;
  45.     Qgs= -ra*V*64512;
  46.     Qrs1= ((Fao*Cpa) + (Fco*Cpc) + (Fto*Cpt))*(T - 310);
  47.     Qrs2= UA*(T - 298);
  48.     Qrs= Qrs1 + Qrs2 + 8.3;
  49.     NCp= (Na*Cpa) + (Nc*Cpc) + (Nt*Cpt) + (Nan*Cpan);
  50.     Caos = 2.22;
  51.     Ccos = 3.52;
  52.     Ctos = 2.92;
  53.     Vos = 0.73;
  54.     vos = 1.5*10**-4;
  55.     UAS = 5;
  56.     Vs = Vos + (vos*t);
  57.     Nas= Cas*Vs;
  58.     Nans= Cans*Vs;
  59.     Faos = Caos*vos;
  60.     Fcos = Ccos*vos;
  61.     Ftos = Ctos*vos;
  62.     ks= (4.01*10**3)*np.exp(-29000/(8.314*Ts))
  63.     ras= -ks*(Cas)**2
  64.     Qgss= -ras*Vs*64512;
  65.     Qrs1s= ((Faos*Cpa) + (Fcos*Cpc) + (Ftos*Cpt))*(Ts - 310);
  66.     Qrs2s= UAS*(Ts - 298);
  67.     Qrss= Qrs1s + Qrs2s + 8.3;
  68.     NCps= (Nas*Cpa) + (Ncs*Cpc) + (Nts*Cpt) + (Nans*Cpan);
  69.     # Differential equations
  70.     dCadt = (vo*10**-5/V)*(Cao - Ca) + ra
  71.     dNcdt = Cco*vo*10**-5
  72.     dNtdt = Cto*vo*10**-5
  73.     dCandt = -ra - (Can*(vo*10**-5/V))
  74.     dTdt=((Qgs-Qrs)/NCp)
  75.     dCasdt = (vos/Vs)*(Caos - Cas) + ras
  76.     dNcsdt = Ccos*vos
  77.     dNtsdt = Ctos*vos
  78.     dCansdt = -ras - (Cans*(vos/V))
  79.     dTsdt=((Qgss-Qrss)/NCps)
  80.     return np.array([dCadt, dNcdt, dNtdt,dCandt, dTdt,dCasdt, dNcsdt, dNtsdt,dCansdt, dTsdt])
  81.  
  82. tspan = np.linspace(0, 3000, 100) # Range for the independent variable
  83. y0 = np.array([Cai,Nci,Nti,0,355,2.12,2.64,2.20,0,355]) # Initial values for the dependent variables
  84.  
  85. #%%
  86. fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
  87. fig.suptitle("""LEP-Synthron Explosion""", fontweight='bold', x = 0.15, y=0.98)
  88. plt.subplots_adjust(left=0.45)
  89. fig.subplots_adjust(wspace=0.25,hspace=0.3)
  90. sol =  odeint(ODEfun, y0, tspan, (Cao,Cco,Cto,Vo,vo,Cai,Nci,Nti,UA))
  91. Ca = sol[:, 0]
  92. Nc= sol[:, 1]
  93. Nt= sol[:, 2]
  94. Can= sol[:, 3]
  95. T= sol[:, 4]
  96. Cas = sol[:, 5]
  97. Ncs= sol[:, 6]
  98. Nts= sol[:, 7]
  99. Cans= sol[:, 8]
  100. Ts= sol[:, 9]
  101. Vos = 0.73;
  102. vos = 1.5*10**-4;
  103. Vs = Vos + (vos*tspan)
  104. V = Vo + (vo*10**-5*tspan);
  105. k= (4.01*10**3)*np.exp(-29000/(8.314*T))
  106. ra= -k*(Ca)**2;
  107. Qgs= -ra*V*64512;
  108. ks= (4.01*10**3)*np.exp(-29000/(8.314*Ts))
  109. ras= -ks*(Cas)**2
  110. Qgss= -ras*Vs*64512
  111. p1,p2= ax2.plot(tspan, Ts,tspan,T)
  112. ax2.legend([r'$T_{Original\hspace{0.5} Recipe}$',r'$T_{Modified \hspace{0.5} Recipe}$'], loc='upper right')
  113. ax2.set_xlabel('time $(sec)$', fontsize='medium')
  114. ax2.set_ylabel(r'$Temp$ (K)', fontsize='medium')
  115. ax2.grid()
  116. ax2.set_ylim(300, 600)
  117. ax2.set_xlim(0, 3000)
  118.  
  119. p3,p4 = ax3.plot(tspan, Cas,tspan,Ca)
  120. ax3.legend([r'$Ca_{\hspace{0.5} Original\hspace{0.5} Recipe}$',r'$Ca_{\hspace{0.5} Modified\hspace{0.5} Recipe}$'], loc='upper right')
  121. ax3.set_xlabel('time $(sec)$', fontsize='medium')
  122. ax3.set_ylabel(r'$Concentration (kmol/m^3)$', fontsize='medium')
  123. ax3.grid()
  124. ax3.set_ylim(0, 0.5)
  125. ax3.set_xlim(0, 3000)
  126.  
  127. p7,p8 = ax4.plot(tspan,Qgss,tspan,Qgs)
  128. ax4.legend([r'$Qg_{\hspace{0.5}  Original\hspace{0.5} Recipe}$',r'$Qg_{\hspace{0.5}  Modified\hspace{0.5} Recipe}$'], loc='upper right')
  129. ax4.set_xlabel('time $(sec)$', fontsize='medium')
  130. ax4.set_ylabel(r'Heat gen $(J/s)$', fontsize='medium')
  131. ax4.grid()
  132. ax4.ticklabel_format(axis='y', style='sci', scilimits=(4,4))
  133. ax4.set_ylim(0, 2e5)
  134. ax4.set_xlim(0, 100)
  135.  
  136. ax1.text(-2,-1.2,'Differential Equations'
  137.          '\n'
  138.          r'$\dfrac{dC_A}{dt} =(v_{0}/V)*(Cao-Ca)+ra$'
  139.          '\n'
  140.          r'$\dfrac{dN_C}{dt} =C_{C0}*v_0$'
  141.          '\n'
  142.          r'$\dfrac{dN_t}{dt} =C_{T0}*v_0$'
  143.          '\n'
  144.          r'$\dfrac{C_{AN}}{dt}=-r_A -(C_{AN}*(v_0/V))$'
  145.          '\n'
  146.           r'$\dfrac{dT}{dt}=(Q_{gs}-Q_{rs})/NCp$'
  147.          '\n\n'
  148.          'Explicit Equations'
  149.          '\n\n'
  150.          r'$V=V_0+(v_0*t)$'
  151.          '\n\n'
  152.          r'$N_{A}=C_A*V$'
  153.         '\n'
  154.          r'$N_{AN}=C_{AN}*V$'
  155.         '\n\n'
  156.          r'$F_{A0}=C_{A0}*v_0$'
  157.          '\n'
  158.          r'$F_{C0}=C_{C0}*v_0$'
  159.          '\n'
  160.          r'$F_{T0}=C_{T0}*v_0$'
  161.          '\n\n'
  162.          r'$k=4.01*10^3*exp(-29000/(8.314*T))$'
  163.          '\n'
  164.          r'$r_{A}=-k*(C_A)^2  $'
  165.          '\n\n'
  166.          r'$Q_{gs}=-r_A*V*64512 $'
  167.          '\n\n'
  168.          r'$Q_{rs1}=(F_{A0}*C_{PA}+F_{C0}*C_{PC}+F_{T0}*C_{PT})*(T-310)$'
  169.          '\n'
  170.          r'$Q_{rs2}=UA*(T-298)$'
  171.          '\n'
  172.          r'$Q_{rs}=Q_{rs1}+Q_{rs2}+8.3$'
  173.          '\n\n'
  174.          r'$N_{CP}=N_{A}*C_{PA}+N_{C}*C_{PC}+N_{T}*C_{PT}+N_{AN}*C_{PAN}$'
  175.           , ha='left', wrap = True, fontsize=13,
  176.         bbox=dict(facecolor='none', edgecolor='black', pad=10.0), fontweight='bold')
  177.  
  178. ax1.axis('off')
  179. axcolor = 'black'
  180. ax_Cao = plt.axes([0.45, 0.82, 0.15, 0.015], facecolor=axcolor)
  181. ax_Vo = plt.axes([0.45, 0.78, 0.15, 0.015], facecolor=axcolor)
  182. ax_vo = plt.axes([0.45, 0.74, 0.15, 0.015], facecolor=axcolor)
  183. ax_Cai = plt.axes([0.45, 0.70, 0.15, 0.015], facecolor=axcolor)
  184. ax_Nci = plt.axes([0.45, 0.66, 0.15, 0.015], facecolor=axcolor)
  185. ax_Nti = plt.axes([0.45, 0.62, 0.15, 0.015], facecolor=axcolor)
  186. ax_UA = plt.axes([0.45, 0.58, 0.15, 0.015], facecolor=axcolor)
  187.  
  188.  
  189. sCao = Slider(ax_Cao, r'$C_{A0}$($\frac{kmol}{m^3}$)', 0.5, 20, valinit=3.99,valfmt='%1.2f')
  190. sVo = Slider(ax_Vo,r'$V_{0}$ ($m^3$)', 0.1,3, valinit= 2.38,valfmt='%1.2f')
  191. svo = Slider(ax_vo,r'$v_0*10^{-5}$ ($\frac{m^3}{s}$)',10**-3, 10**3, valinit=6,valfmt='%1.2f')
  192. sCai = Slider(ax_Cai,r'$C_{A}(0)$($\frac{kmol}{m^3}$)', 0.5, 10, valinit= 2.38,valfmt='%1.2f')
  193. sNci = Slider(ax_Nci,r'$N_{C}(0)$ ($kmol$)', 0.5, 20, valinit= 7.76,valfmt='%1.2f')
  194. sNti = Slider(ax_Nti,r'$N_{T}(0)$ ($kmol$)', 0.5, 20, valinit= 7.09,valfmt='%1.2f')
  195. sUA = Slider(ax_UA,r'$UA$ ($\frac{J}{s.K}$)', 0,15, valinit= 5,valfmt='%1.1f')
  196.  
  197. def update_plot2(val):
  198.     Cao = sCao.val
  199.     Vo =sVo.val
  200.     vo = svo.val
  201.     Cai =sCai.val
  202.     Nci = sNci.val
  203.     Nti =sNti.val
  204.     UA =sUA.val
  205.     y0 = np.array([Cai,Nci,Nti,0,355,2.12,2.64,2.20,0,355])
  206.     sol = odeint(ODEfun, y0, tspan, (Cao,Cco,Cto,Vo,vo,Cai,Nci,Nti,UA))
  207.     Ca = sol[:, 0]
  208.     T= sol[:, 4]
  209.     Cas = sol[:, 5]
  210.     Ts= sol[:, 9]
  211.     V = Vo + (vo*10**-5*tspan);
  212.     k= (4.01*10**3)*np.exp(-29000/(8.314*T))
  213.     ra= -k*(Ca)**2;
  214.     Qgs= -ra*V*64512;
  215.  
  216.     p1.set_ydata(Ts)
  217.     p2.set_ydata(T)
  218.     p3.set_ydata(Cas)
  219.     p4.set_ydata(Ca)
  220.     p7.set_ydata(Qgss)
  221.     p8.set_ydata(Qgs)
  222.     fig.canvas.draw_idle()
  223.  
  224. sCao.on_changed(update_plot2)
  225. sVo.on_changed(update_plot2)
  226. svo.on_changed(update_plot2)
  227. sCai.on_changed(update_plot2)
  228. sNci.on_changed(update_plot2)
  229. sNti.on_changed(update_plot2)
  230. sUA.on_changed(update_plot2)
  231.  
  232. resetax = plt.axes([0.48, 0.87, 0.09, 0.04])
  233. button = Button(resetax, 'Reset variables', color='cornflowerblue', hovercolor='0.975')
  234.  
  235. def reset(event):
  236.     sCao.reset()
  237.     sVo.reset()
  238.     svo.reset()
  239.     sCai.reset()
  240.     sNci.reset()
  241.     sNti.reset()
  242.     sUA.reset()
  243.  
  244. button.on_clicked(reset)
Raw Paste

Login or Register to edit or fork this paste. It's free.