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))
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)
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,
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

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