Non-Isothermal CSTR Simulation
By Taylor Geisler
The Problem
In this example we will investigate a non-isothermal stirred tank reactor with a single, nonreversible reaction.
\(A+B\rightarrow C+D\)
The concentration over time of each species can be described by an ordinary differential equation. The ODE consists of terms accounting for flow in, flow out, generation, and consumption.
\(\frac{dN_{i}}{dt}=flow\: in-flow\: out+generation-consumption\)
Substituting the rate of reaction for the the generation and consumption terms we get
\(\frac{dN_{i}}{dt}=F_{i\, in}-F_{i\, out}+r_{i}V\)
Dividing by the volume of the reactor V we get the differential equation in terms of concentration.
\(\frac{dC_{i}}{dt}=\frac{\nu_{in}}{V}C_{i\, in}-\frac{\nu_{out}}{V}C_{i\, out}+r_{i}\)
Using this ODE we can solve for the reactor contents over time. The equations for each compenent are written here.
\(\frac{dC_{A}}{dt}=\frac{\nu_{in}}{V}C_{A\, in}-\frac{\nu_{out}}{V}C_{A\, out}+r_{A}\)
\(\frac{dC_{B}}{dt}=\frac{\nu_{in}}{V}C_{B\, in}-\frac{\nu_{out}}{V}C_{B\, out}+r_{B}\)
\(\frac{dC_{C}}{dt}=\frac{\nu_{in}}{V}C_{C\, in}-\frac{\nu_{out}}{V}C_{C\, out}+r_{C}\)
\(\frac{dC_{D}}{dt}=\frac{\nu_{in}}{V}C_{D\, in}-\frac{\nu_{out}}{V}C_{D\, out}+r_{D}\)
The temperature is also described by a differential equation.
\(\frac{dT}{dt}=\frac{Q-\sum_{i=1}^{n}F_{i0}(C_{Pi}(T-T_{i0}))-\Delta H_{RX}(-r_{A}V)}{\sum_{i=1}^{n}N_{i}C_{Pi}}\)
We will use the common irreversible rate equation
\(r_{A}=kC_{A}C_{B}\)
With k as a function of temperature.
\(k=k_{1}(T_{1})exp\left[\frac{E}{R}(\frac{1}{T_{1}}-\frac{1}{T})\right]\)
And
\(r_{A}=r_{B}=-r_{C}=-r_{D}\)
We must know the initial concentrations and temperature in the reactor in order to solve the differential equations. Once we know these we have everything we need to solve for the concentration over time.
The Solver
The Backward Euler method was used to solve the ordinary differential equations. This is an implicit time integrator that uses the derivative at time step j+1 to determine the change in concentration from time j to j+1.
\(C_{i\, j+1}=C_{i\, j}+\Delta t\frac{dC_{i\, j+1}}{dt}\)
\(T_{j+1}=T_{j}+\Delta t\frac{dT_{j+1}}{dt}\)
Here is this equation written in residual form for each species. Notice that the differential expressions from above have been substituted for the differential.
\(f_{1}=-C_{A\, j+1}+C_{A\, j}+\Delta t\left(\frac{\nu_{in}}{V}C_{A\, in}-\frac{\nu_{out}}{V}C_{A\, j+1}+r_{A\, j+1}\right)\)
\(f_{2}=-C_{B\, j+1}+C_{B\, j}+\Delta t\left(\frac{\nu_{in}}{V}C_{B\, in}-\frac{\nu_{out}}{V}C_{B\, j+1}+r_{B\, j+1}\right)\)
\(f_{3}=-C_{C\, j+1}+C_{C\, j}+\Delta t\left(\frac{\nu_{in}}{V}C_{C\, in}-\frac{\nu_{out}}{V}C_{C\, j+1}+r_{C\, j+1}\right)\)
\(f_{4}=-C_{D\, j+1}+C_{D\, j}+\Delta t\left(\frac{\nu_{in}}{V}C_{D\, in}-\frac{\nu_{out}}{V}C_{D\, j+1}+r_{D\, j+1}\right)\)
\(f_{5}=-T_{j+1}+T_{j}+\Delta t\frac{Q-\sum_{i=1}^{n}F_{i0}(C_{Pi}(T-T_{i0}))-\Delta H_{RX}(-r_{A}V)}{\sum_{i=1}^{n}N_{i}C_{Pi}}\)
We cannot solve these equations algebraically. Note that the r term is a nonlinear function of the product of \( C_{A} \) and \( C_{B} \).
To solve this nonlinear system of equations, Newton’s Method was employed at each time step. Newton’s method requires iteratively solving the equation
\(J\cdot\Delta x=-RHS\)
The terms J and RHS are constructed using the nonlinear equations that we are trying to solve.
\(
J=\begin{bmatrix}\frac{df_{1}}{dC_{A}} & \frac{df_{1}}{dC_{B}} & \frac{df_{1}}{dC_{C}} & \frac{df_{1}}{dC_{D}} & \frac{df_{1}}{dT}\\
\frac{df_{2}}{dC_{A}} & \frac{df_{2}}{dC_{B}} & \frac{df_{2}}{dC_{C}} & \frac{df_{2}}{dC_{D}} & \frac{df_{2}}{dT}\\
\frac{df_{3}}{dC_{A}} & \frac{df_{3}}{dC_{B}} & \frac{df_{3}}{dC_{C}} & \frac{df_{3}}{dC_{D}} & \frac{df_{3}}{dT}\\
\frac{df_{4}}{dC_{A}} & \frac{df_{4}}{dC_{B}} & \frac{df_{4}}{dC_{C}} & \frac{df_{4}}{dC_{D}} & \frac{df_{4}}{dT}\\
\frac{df_{5}}{dC_{A}} & \frac{df_{5}}{dC_{B}} & \frac{df_{5}}{dC_{C}} & \frac{df_{5}}{dC_{D}} & \frac{df_{5}}{dT}
\end{bmatrix}
\)
\(
RHS=\begin{bmatrix}f_{1}(C_{A},C_{B},C_{C},C_{D,}T)\\
f_{2}(C_{A},C_{B},C_{C},C_{D,}T)\\
f_{3}(C_{A},C_{B},C_{C},C_{D},T)\\
f_{4}(C_{A},C_{B},C_{C},C_{D},T)\\
f_{5}(C_{A},C_{B},C_{C},C_{D},T)
\end{bmatrix}
\)
\(
\Delta x=\begin{bmatrix}\Delta C_{A}\\
\Delta C_{B}\\
\Delta C_{C}\\
\Delta C_{D}\\
\Delta T
\end{bmatrix}
\)
First, a guess of \( C_{A\, j+1}\), \( C_{B\, j+1}\), \( C_{C\, j+1}\),\( C_{D\, j+1} \) and T are required. The jacobian matrix J and the RHS vector are calculated. The derivatives in the jacobian matrix were calculated using the forward difference approximation.
Now that we have all the necessary terms, we can solve for the vector \( \Delta x\). This is a vector of updates to our initial guesses. We add these updates to our original guesses and reconstruct the jacobian and RHS terms with the new, better values. This process is repeated until the updates become small. This means that we have converged to the solution, within a given tolerance.
By repeating this process at each time step, we solve for \( C_{A}\), \( C_{B}\), \( C_{C}\), \( C_{D} \), and T over time. The solution is plotted above.
The Simulation
Experiment with different values in the above simulation to see how the reaction will change. The simulation shows both inital reactor dynamics as well as the steady state values that are eventually reached. Try different values of the inital temperature to change the reaction rate. Experiment with changing the rate constant of the reaction to produce a higher or lower concentration of products at steady state. Adjust the flow rate and examine the change in the inital dynamics of the system. Examine how the steady state values react when you change the inlet concentration of one of the products or reactants. The heat into the system can also be adjusted.