# This code is an extension of http://scpyce.org/12/importnumpyasnpfromscipyimportintegratefrommatplotlib.pylabimport*deftank(t,y):""" Dynamic balance for a CSTR C_A = y[0] = the concentration of A in the tank, [mol/L] T = y[1] = the tank temperature, [K] Returns dy/dt = [F/V*(C_{A,in} - C_A) - k*C_A^2 ] [F/V*(T_in - T) - k*C_A^2*HR/(rho*Cp) ] """F=20.1# L/minCA_in=2.5# mol/LV=100.0# Lk0=0.15# L/(mol.min)Ea=5000# J/molR=8.314# J/(mol.K)Hr=-590# J/molT_in=288# Krho=1.050# kg/L# Assign some variables for convenience of notationCA=y[0]T=y[1]# Algebraic equationsk=k0*np.exp(-Ea/(R*T))# L/(mol.min)Cp=4.184-0.002*(T-273)# J/(kg.K)# Output from ODE function must be a COLUMN vector, with n rowsn=len(y)# 2: implies we have two ODEsdydt=np.zeros((n,1))dydt[0]=F/V*(CA_in-CA)-k*CA**2dydt[1]=F/V*(T_in-T)-(Hr*k*CA**2)/(rho*Cp)returndydt# The "driver" that will integrate the ODE(s):# -----------# Start by specifying the integrator:# use ``vode`` with "backward differentiation formula"r=integrate.ode(tank).set_integrator('vode',method='bdf')# Set the time ranget_start=0.0t_final=45.0delta_t=0.1# Number of time steps: 1 extra for initial conditionnum_steps=np.floor((t_final-t_start)/delta_t)+1# Set initial condition(s): for integrating variable and time!CA_t_zero=0.5T_t_zero=295.0r.set_initial_value([CA_t_zero,T_t_zero],t_start)# Additional Python step: create vectors to store trajectoriest=np.zeros((num_steps,1))CA=np.zeros((num_steps,1))temp=np.zeros((num_steps,1))t[0]=t_startCA[0]=CA_t_zerotemp[0]=T_t_zero# Integrate the ODE(s) across each delta_t timestepk=1whiler.successful()andk<num_steps:r.integrate(r.t+delta_t)# Store the results to plot latert[k]=r.tCA[k]=r.y[0]temp[k]=r.y[1]k+=1# All done! Plot the trajectories in two separate plots:fig=figure()ax1=subplot(211)ax1.plot(t,CA)ax1.set_xlim(t_start,t_final)ax1.set_xlabel('Time [minutes]')ax1.set_ylabel('Concentration [mol/L]')ax1.grid('on')ax2=plt.subplot(212)ax2.plot(t,temp,'r')ax2.set_xlim(t_start,t_final)ax2.set_xlabel('Time [minutes]')ax2.set_ylabel('Temperaturere [K]')ax2.grid('on')fig.savefig('coupled-ode-plot.png')

Software license:Creative Commons Zero. No rights reserved.
Users have permission to do anything with the code and other material on this page. (More details)

More information:

We extend the example from http://scpyce.org/12/ here to integrate 2 coupled ODEs and include several algebraic equations.

Consider a chemical reactor with a second order reaction \({\sf A} \rightarrow {\sf B}\). Our aim is calculate the concentration, \(C_{\sf A}(t)\), and tank temperature, \(T(t)\), as functions of time using an ODE integrator.

We know the reaction rate is \(r = k C_{\sf A}^2\), and is an algebraic function of temperature: \(k = 0.15 e^{- E_a/(RT)}\) L/(mol.min). Furthermore, the reaction is known to release heat, with \(\Delta H_r = -590\) J/mol. The time-dependent mass and energy balance for this system can be derived:

In the code we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 45.0\) minutes
and plot the time-varying trajectories of temperature and concentration in the tank using matplotlib.