Appendix I — Axial Dispersion Reactor Models
This appendix starts with the mole and energy balances for an ideal PFR. It adds dispersion terms to allow for mixing of heat and mass in the axial direction, leading to one form of the transient axial dispersion reactor design equations. In anticipation of solving the steady-state versions of those equations numerically, it then shows how the steady state axial dispersion mole and energy balances can be re-written as an equivalent set of first-order differential equations. Conventionally the axial dispersion model is written using the superficial velocity of the fluid in the \(z\)-direction, \(u_z\), but for consistency with the other chapters in Reaction Engineering Basics it is written here in terms of the volumetric flow rate.
I.1 Axial Dispersion Reactor Design Equations
The transient mole and energy balances for a PFR are derived in Appendix H. They serve as the starting point for the derivation of the axial dispersion reactor design equations. The assumption of plug flow is retained in the axial dispersion model. The only modification to the PFR mole and energy balances is addition of the term, \(-D_{ax} \frac{\pi D^2}{4} \frac{d^2C_i}{dz^2}\), to the ideal PFR mole balance, and the term, \(-\lambda_{ax} \frac{\pi D^2}{4} \frac{d^2T}{dz^2}\) to the ideal PFR reacting fluid energy balance.
I.1.1 The Transient Axial Dispersion Mole Balance
The generation of the axial dispersion mole balance begins with ideal PFR mole balance.
\[ \frac{\pi D^2}{4\dot V} \frac{\partial\dot n_i}{\partial t} - \frac{\pi D^2\dot n_i}{4\dot V^2} \frac{\partial \dot V}{\partial t} + \frac{\partial \dot n_i}{\partial z} = \frac{\pi D^2}{4}\sum_j \nu_{i,j}r_j \]
The dispersion term is added, noting that because \(C_i\) varies with both \(t\) and \(z\), it uses a partial derivative.
\[ \frac{\pi D^2}{4\dot V} \frac{\partial\dot n_i}{\partial t} - \frac{\pi D^2\dot n_i}{4\dot V^2} \frac{\partial \dot V}{\partial t} - D_{ax} \frac{\pi D^2}{4} \frac{\partial ^2C_i}{\partial z^2} + \frac{\partial \dot n_i}{\partial z} = \frac{\pi D^2}{4}\sum_j \nu_{i,j}r_j \]
This presents a dilemma because the concentration, \(C_i\), is the dependent variable in the axial dispersion term, but the molar flow rate, \(\dot{n}_i\), is the dependent variable in the ideal PFR mole balance. Consequently, either the dispersion term must be re-written in terms of the molar flow rate, or, the molar flow rate of the PFR mole balance must be re-written in terms of the concentration. Practically, the axial dispersion mole balance is much more tractable when it is written using the concentration as the dependent variable, so the terms including the molar flow rate need to be re-expressed in terms of the concentration. To do so, both sides of the equation first are multiplied by by \(\frac{4}{\pi D^2}\).
\[ \frac{1}{\dot V} \frac{\partial\dot n_i}{\partial t} - \frac{\dot n_i}{\dot V^2} \frac{\partial \dot V}{\partial t} - D_{ax} \frac{\partial ^2C_i}{\partial z^2} + \frac{4}{\pi D^2} \frac{\partial \dot n_i}{\partial z} = \sum_j \nu_{i,j}r_j \]
Noting that \(C_i = \frac{\dot{n}_i}{\dot{V}}\) so \(\frac{\partial C_i}{\partial t} = \frac{1}{\dot V} \frac{\partial\dot n_i}{\partial t} - \frac{\dot n_i}{\dot V^2} \frac{\partial \dot V}{\partial t}\) facilitates conversion of the dependent variable in the time derivatives.
\[ \frac{\partial C_i}{\partial t} - D_{ax} \frac{\partial ^2C_i}{\partial z^2} + \frac{4}{\pi D^2} \frac{\partial \dot n_i}{\partial z} = \sum_j \nu_{i,j}r_j \]
Similarly noting that \(\dot{n}_i\) = \(\dot{V} C_i\) eliminates the molar flow rate as the dependent varibles in the \(z\) derivative.
\[ \frac{\partial C_i}{\partial t} - D_{ax} \frac{\partial ^2C_i}{\partial z^2} + \frac{4}{\pi D^2} \frac{\partial }{\partial z}\left( \dot{V} C_i \right) = \sum_j \nu_{i,j}r_j \]
Applying the chain rule for the derivative of a product, \(\frac{\partial }{\partial z}\left( \dot{V} C_i \right) = \dot{V} \frac{\partial C_i}{\partial z} + C_i \frac{\partial \dot{V}}{\partial z}\) then yields the desired form of the transient axial dispersion mole balance.
\[ \frac{\partial C_i}{\partial t} - D_{ax}\frac{\partial^2 C_i}{\partial z^2} + \frac{4 \dot{V}}{\pi D^2}\frac{\partial C_i}{\partial z} + \frac{4 C_i}{\pi D^2}\frac{\partial \dot{V}}{\partial z} = \sum_j \nu_{i,j}r_j \]
I.1.2 The Transient Axial Dispersion Energy Balance
The generation of the energy balance begins with the ideal PFR reacting fluid energy balance.
\[ \begin{align} \frac{\pi D^2}{4 \dot{V}} \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)& \frac{\partial T}{\partial t} - \frac{\pi D^2}{4} \frac{\partial P}{\partial t} + \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)\frac{\partial T}{\partial z} \\&= \pi DU\left( T_{ex} - T \right) - \frac{\pi D^2}{4} \sum_j r_j\Delta H_j \end{align} \]
The dispersion term is added, again noting that because \(T\) varies with both \(t\) and \(z\), a partial derivative is used.
\[ \begin{align} \frac{\pi D^2}{4 \dot{V}} \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)& \frac{\partial T}{\partial t} - \frac{\pi D^2}{4} \frac{\partial P}{\partial t} -\lambda_{ax} \frac{\pi D^2}{4} \frac{d^2T}{dz^2} + \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)\frac{\partial T}{\partial z} \\&= \pi DU\left( T_{ex} - T \right) - \frac{\pi D^2}{4} \sum_j r_j\Delta H_j \end{align} \]
Both sides of the equation can then be multiplied by \(\frac{4}{\pi D^2}\).
\[ \begin{align} \frac{1}{\dot{V}} \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)& \frac{\partial T}{\partial t} - \frac{\partial P}{\partial t} - \lambda_{ax} \frac{d^2T}{dz^2} + \frac{4}{\pi D^2} \left( \sum_i \dot{n}_i\hat{C}_{p,i} \right)\frac{\partial T}{\partial z} \\&= \frac{4U}{D} \left( T_{ex} - T \right) - \sum_j r_j\Delta H_j \end{align} \]
Noting that \(\dot{n}_i\) = \(\dot{V} C_i\) allows elimination of the molar flow rate, and after doing so, \(\dot{V}\) can be factored out of the summations yielding the desired form of the axial dispersion energy balance.
\[ \begin{align} \left( \sum_i C_i\hat{C}_{p,i} \right)& \frac{\partial T}{\partial t} - \frac{\partial P}{\partial t} - \lambda_{ax} \frac{d^2T}{dz^2} + \frac{4\dot{V}}{\pi D^2} \left( \sum_i C_i\hat{C}_{p,i} \right)\frac{\partial T}{\partial z} \\&= \frac{4U}{D} \left( T_{ex} - T \right) - \sum_j r_j\Delta H_j \end{align} \]
I.1.3 The Steady-State Axial Dispersion Reactor Design Equations.
In this book transient, non-ideal, tubular reactors will not be analyzed. The transient design equations are presented here for the sake of completeness. When the reactor operates at steady-state, all of the time derivatives are equal to zero, and the spatial derivatives become ordinary derivatves.
\[ - D_{ax}\frac{d^2 C_i}{d z^2} + \frac{4 \dot{V}}{\pi D^2}\frac{d C_i}{d z} + \frac{4 C_i}{\pi D^2}\frac{d \dot{V}}{d z} = \sum_j \nu_{i,j}r_j \]
\[ - \lambda_{ax}\frac{d^2 T}{d z^2} + \frac{4 \dot{V}}{\pi D^2} \left(\sum_i C_i \hat{C}_{p,i}\right) \frac{d T}{d z} = \frac{4U}{D}\left(T_{ex} - T\right) - \sum_j r_j \Delta H_j \]
The mole and energy balances are second-order ordinary differential equations, and the number of dependent variables is one greater than the number of equations. This is due to the presence of the derivative of the volumetric flow rate with respect to \(z\). As a consequence, either one dependent variable must be eliminated, or another ODE must be added. Noting that mass is neither created nor destroyed in the reactor a steady-state total mass balance can be used to provide the added ODE.
\[ \frac{d \dot{m}}{dz} = 0. \]
The mass flow rate is equal to the product of the volumetric flow rate and the fluid density.
\[ \frac{d}{dz}\left( \rho \dot{V}\right) = \rho \frac{d\dot{V}}{dz} + \dot{V} \frac{d\rho}{dz} = 0. \]
Rearranging yields an expression for the derivative of the volumetric flow rate in terms of the fluid density.
\[ \frac{d\dot{V}}{dz} = \frac{-1}{\rho}\frac{d\rho}{dz} \]
In general, the density is not a constant, but it can be expressed in terms of the concentrations and molecular weights of the reagents.
\[ \rho = \sum_i \left(M_i C_i\right) \qquad \Rightarrow \qquad \frac{d \rho}{dz} = \sum_i \left(M_i \frac{dC_i}{dz}\right) \]
Substitution yields the desired ODE.
\[ \frac{d\dot{V}}{dz} = \frac{-\sum_i \left(M_i \frac{dC_i}{dz}\right)}{\sum_i \left(M_i C_i\right)} \]
I.1.4 Steady-State Axial Dispersion Model Boundary Conditions
The axial dispersion design conditions include second derivatives of each concentration and of temperature and the first derivative of the volumetric flow rate. Solving them then requires two boundary conditions for each concentration and for the temperature and one boundary condition for the volumetric flow rate. There are a few different ways to write the boundary conditions. The Danckwerts boundary conditions will be presented here.
It is important to note that the concentrations and temperature at the inlet to the reactor (i. e. at \(z=0\)) will not equal the feed concentrations and temperature, as indicated in Figure I.1. That is, due axial dispersion, there is mixing at the reactor inlet between the feed flow and the fluid within the reactor.
The Danckwerts boundary conditions account for that mixing. (Here it has been assumed that the heat capacities are constant.)
\[ \dot{V}_{feed}C_{i,feed} - D_{ax} \frac{\pi D^2}{4} \frac{\partial C_i}{\partial z}\Bigg\vert_{z=0} = \dot{V}_0C_{i,0} \]
\[ \dot{V}_{feed}\sum_i\left(C_{i,feed} \hat{C}_{p,i} \right) T_{feed} - \lambda_{ax} \frac{\pi D^2}{4} \frac{\partial T}{\partial z}\Bigg\vert_{z=0} = \dot{V}_{0} \sum_i\left(C_{i,0} \hat{C}_{p,i} \right) T_0 \]
The volumetric flow rate at the inlet can be found by noting that the mass flow rate is constant.
\[ \dot{m}_{0} = \dot{m}_{feed} \qquad \Rightarrow \qquad \rho_0 \dot{V}_0 = \rho_{feed} \dot{V}_{feed} \qquad \Rightarrow \qquad \dot{V}_0 = \frac{\rho_{feed}}{\rho_0}\dot{V}_{feed} \]
\[ \dot{V}_0 = \frac{\sum_i \left(M_iC_{i,feed}\right)}{\sum_i\left(M_iC_{i,0}\right)}\dot{V}_{feed} = \frac{\dot{m}_{feed}}{\sum_i\left(M_iC_{i,0}\right)} \]
A second boundary condition is needed for each \(C_i\) and for \(T\). Because the reactor ends at \(z=L\), it can be assumed that the concentrations and temperature will stop changing at that point. Doing so yields the additional boundary conditions.
\[ \frac{\partial C_i}{\partial z}\Bigg\vert_{z=L} = 0 \]
\[ \frac{\partial T}{\partial z}\Bigg\vert_{z=L} = 0 \]
I.2 Equivalent First-Order, Steady-State Design Equations and Boundary Conditions
Recognizing that the design equations will be solved numerically, it is useful to rewrite them as an equivalent set of first-order ODEs. This can be accomplished by defining new variables, \(\omega_i\) and \(\omega_T\).
\[ \frac{dC_i}{dz} = \omega_i \]
\[ \frac{dT}{dz} = \omega_T \]
Noting that then \(\frac{d^2C_i}{dz^2} = \frac{d\omega_i}{dz}\) and \(\frac{d^2T}{dz^2} = \frac{d\omega_T}{dz}\), substitution into the second-order design equations converts them to first-order ODEs.
\[ - D_{ax}\frac{d \omega_i}{d z} + \frac{4 C_i}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_i + \sum_j \nu_{i,j}r_j \]
\[ - \lambda_{ax}\frac{d \omega_T}{d z} = - \frac{4 \dot{V}}{\pi D^2} \left(\sum_i C_i \hat{C}_{p,i}\right) \omega_T + \frac{4U}{D}\left(T_{ex} - T\right) - \sum_j r_j \Delta H_j \]
The boundary conditions must also be expressed in terms of the new variables. Substitution in the Danckwerts boundary conditions and rearrangement yields their equivalent first-order forms.
\[ C_{i,0} = \frac{\dot{V}_{feed}C_{i,feed} - D_{ax} \frac{\pi D^2}{4} \omega_{i,0}}{\dot{V}_0} \]
\[ T_0 = \frac{\dot{V}_{feed}\sum_i\left(C_{i,feed} \hat{C}_{p,i} \right) T_{feed} - \lambda_{ax} \frac{\pi D^2}{4} \omega_{T,0}}{\dot{V}_{0} \sum_i\left(C_{i,0} \hat{C}_{p,i} \right) } \]
Similarly, substitution in the boundary conditions at \(z=L\) yields the equivalent boundary conditions for the first-order design equations.
\[ \omega_i\Big\vert_{z=L} = \omega_{i,1} = 0 \]
\[ \omega_T\Big\vert_{z=L} = \omega_{T,1} = 0 \]
I.3 Axial Dispersion Design Equations as CVODEs
Some of the boundary conditions are specified at the reactor inlet and some of them at the reactor outlet. In particular, the values of each \(\omega_i\) and of \(\omega_T\) are not known at the reactor inlet, and they are needed to calculate the values of each \(C_i\) and of \(T\) at the inlet. At the same time, the values of each \(\omega_i\) and of \(\omega_T\) are known at the reactor outlet. This makes it possible to solve the axial dispersion reactor design equations using an IVODE solver. That is, the reactor design equations are CVODEs (see Appendix M).
To solve the axial dispersion design equations using an IVODE solver, the unknown initial values of each \(\omega_i\), \(\omega_T\), and \(\dot{V}\) are considered coupled values. The known values of each \(\omega_i\) and \(\omega_T\) at the outlet are used to write implicit equations for the inlet values together with he explicit equation for \(\dot{V}_0\) from above.
\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{i,1} = 0 \]
\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{T,1} = 0 \]
Effectively, these equations mean find the values of each \(\omega_i\) and \(\omega_T\) at the inlet such that their values equal zero at the outlet.
Then, when it is necessary to solve the the axial dispersion reactor design equations, these implicit equations can first solved using a ATE solver. Doing so yields the values of each \(\omega_i\) and \(\omega_T\) at the inlet with which the inlet values of each \(C_i\) and of \(T\) can be calculated. Then, knowing the inlet values of every dependent variable, the axial dispersion reactor design equations can be solved using an IVODE solver.
I.4 Summary
The axial dispersion mole and energy balances can be written as the following set of first-order CVODEs. If there is pressure drop in the reactor the PFR momentum balance should be added to the design equations. If the design equations include an energy balance on a heat exchange fluid, it will be an ATE that can be solved simultaneously with the coupled value equations, if necessary.
\[ - D_{ax}\frac{d \omega_i}{d z} + \frac{4 C_i}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_i + \sum_j \nu_{i,j}r_j \]
\[ - \lambda_{ax}\frac{d \omega_T}{d z} = - \frac{4 \dot{V}}{\pi D^2} \left(\sum_i C_i \hat{C}_{p,i}\right) \omega_T + \frac{4U}{D}\left(T_{ex} - T\right) - \sum_j r_j \Delta H_j \]
\[ \frac{d\dot{V}}{dz} = \frac{-\sum_i \left(M_i \frac{dC_i}{dz}\right)}{\sum_i \left(M_i C_i\right)} \]
\[ \frac{dC_i}{dz} = \omega_i \]
\[ \frac{dT}{dz} = \omega_T \]
The initial values of each \(\omega_i\) and of \(\omega_T\) are coupled values that can be found by solving the following implicit coupled value equations using an ATE solver.
\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{i,1} = 0 \]
\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{T,1} = 0 \]
The remaining initial values then can be calculated as follows.
\[ C_{i,0} = \frac{\dot{V}_{feed}C_{i,feed} - D_{ax} \frac{\pi D^2}{4} \omega_{i,0}}{\dot{V}_0} \]
\[ T_0 = \frac{\dot{V}_{feed}\sum_i\left(C_{i,feed} \hat{C}_{p,i} \right) T_{feed} - \lambda_{ax} \frac{\pi D^2}{4} \omega_{T,0}}{\dot{V}_{0} \sum_i\left(C_{i,0} \hat{C}_{p,i} \right) } \]
\[ \dot{V}_0 = \frac{\sum_i \left(M_iC_{i,feed}\right)}{\sum_i\left(M_iC_{i,0}\right)}\dot{V}_{feed} = \frac{\dot{m}_{feed}}{\sum_i\left(M_iC_{i,0}\right)} \]
Then, knowing the inlet values of every dependent variable, the axial dispersion design equations can be solved using an IVODE solver. If there is pressure drop in the reactor, the appropriate PFR momentum balance should be added to the design equations and the inlet pressure should be added to the initial values. If an exchange fluid energy balance is included in the design equations, it can be added to the coupled value equations and solved to find \(T_{ex}\).
Finally, if the reacting fluid is an ideal, incompressible liquid mixture, the derivative of the volumetric flow rate will equal zero and an initial value for the volumetric flow rate will not be needed.
I.5 Symbols Used in Appendix I
| Symbol | Meaning |
|---|---|
| \(t\) | Time. |