22 The Steady-State Axial Dispersion Model

The ideal PFR model assumes there is no mixing in the axial direction. This chapter presents the steady-state version of a model wherein terms have been added to the PFR mole and energy balances to allow mixing in the axial direction. It is called the axial dispersion model. The design equations for the axial dispersion model are derived from the PFR design equations in Appendix I. It should be noted that the axial dispersion model is usually written using the linear velocity of the fluid in the \(z\) direction, \(u_z\), but for consistenccy with the other chapters of Reaction Engineering Basics, it is written in terms of the volumetric flow rate, \(\dot{V}\), in this chapter.

22.1 Axial Dispersion

A tubular or packed bed reactor only approaches plug flow at large Reynolds numbers where the flow is highly turbulent. As a consequence, there are small eddies and currents that flow in the axial direction at velocities above and below the mean fluid velocity. These small eddies cause some mixing in the axial direction. To account for this, the term, \(-D_{ax} \frac{\pi D^2}{4} \frac{d^2C_i}{dz^2}\), can be added to the ideal PFR mole balance, and the term, \(-\lambda_{ax} \frac{\pi D^2}{4} \frac{d^2T}{dz^2}\), can be added to the PFR energy balance. To those familiar with heat and mass transfer, these terms look terms for mass diffusion and thermal conduction in the \(z\) direction, but they are not.

The coefficient, \(D_{ax}\), is not a diffusion coefficient. It is called a dispersion coefficient, and it is typically obtained from correlations (in dimensionless form) as an axial or longitudinal Peclet number, Equation 22.1. Unlike diffusion coefficients, for a given flow stream, there is a single axial dispersion coefficient that applies for every reagent in the system. Similarly, the coefficient, \(\lambda_{ax}\), is not a thermal condictivity, it is a thermal dispersion coefficient. It can be found in correlations as a dimensionless longitudinal Bodenstein number, Equation 22.2. As noted, the axial dispersion and thermal dispersion coefficients are meant to account for axial mixing due to the turbulence of the flow. Typically, axial mixing of this type is only significant in extremely short tubular reactors.

\[ Pe_{ax} = \frac{4\dot{V}L}{\pi D^2 D_{ax}} \tag{22.1}\]

\[ Bo_{ax} = \frac{4\rho \tilde{C}_p \dot{V}L}{\pi D^2 \lambda_{ax}} \tag{22.2}\]

So, when used to account for axial mixing due to turbulence, the axial dispersion terms are only necessary if the reactor is very short. However, if the flow in a tubular reactor, and particularly in a packed bed reactor, is irregular, the axial dispersion model can be used as an empirical reactor model. In this situation, \(D_{ax}\) and \(\lambda_{ax}\) are used as adjustable parameters whose values are chosen to make the reactor model as accurate as possible. Generally when \(Pe_{ax}\) is less than 10 and \(Bo_{ax}\) is less than 1, backmixing is significant, and when \(Pe_{ax}\) and \(Bo_{ax}\) are greater than 100, backmixing is minimal. Thus, in an empirical reactor model, \(D_{ax}\) and \(\lambda_{ax}\) become parameters that are adjusted to tune the amount of backmixing.

22.2 Axial Dispersion Reactor Design Equations

As noted, the axial dispersion model is simply the ideal PFR model with dispersion terms added to the mole and energy balances. Importantly, plug flow is still assumed. In Appendix I the axial dispersion mole and energy balances and the associated initial values and boundary conditions are derived, the steady-state forms are generated, and then they are converted to an equivalent set of first-order CVODEs. Familiarity with the first-order design equations and the variables (\(\omega_i\) and \(\omega_T\)) introduced with them is assumed here.

Figure I.1 presents a schematic representation of an adiabatic reactor that can be modeled using the axial dispersion model. It is important to note that the concentrations and temperature at the inlet to the reactor will not equal the feed concentrations and temperature. This is a result of axial mixing occurring at the reactor inlet thereby causing some of the feed to mix with fluid already in the reactor.

Figure 22.1: Due to mixing in the axial direction, the concentrations and temperature at the inlet are not equal to the feed concentrations and temperature.

It becomes important to differentiate between the values of variables in the feed, at the reactor inlet, and at the reactor outlet. Here, if the location is not identified using a vertical bar and subscript, e. g. \(C_A \big\vert_{z=0}\) for the concentration of reagent A at the reactor inlet, an additional subscript is added to the variable symbol. Thus, \(C_{A,feed}\) denotes the concentration of A in the feed, \(C_A \big\vert_{z=0}\) or \(C_{A,0}\) denotes the concentration of A at the reactor inlet, and \(C_A \big\vert_{z=L}\) or \(C_{A,1}\) denotes the concentration at the reactor outlet.

The mole and energy balances for the axial dispersion model are shown in Equations 22.3 through 22.7. Equations 22.3 and 22.4 are the PFR design equations with the dispersion terms added. Equation 22.5 is the total mass balance, re-written in terms of the volumetric flow rate. If the reacting fluid is an ideal, incompressible liquid mixture, the volumetric flow rate is constant and equal to \(\dot{V}_{feed}\), the derivative of the volumetric flow rate in Equation 22.3 equals zero, and Equation 22.5 is not needed. Equations 22.6 and 22.7 are essentially defining equations for \(\omega_i\) and \(\omega_T\). If there is pressure drop in the reactor, a momentum balance must be added to the design equations. If an external heat exchange fluid transfers sensible heat to or from the reactor, an energy balance on the exchange fluid must be added, but it will be an ATE, not an ODE.

\[ - 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 \tag{22.3}\]

\[ - \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 \tag{22.4}\]

\[ \frac{d\dot{V}}{dz} = \frac{-\sum_i \left(M_i \frac{dC_i}{dz}\right)}{\sum_i \left(M_i C_i\right)} \tag{22.5}\]

\[ \frac{dC_i}{dz} = \omega_i \tag{22.6}\]

\[ \frac{dT}{dz} = \omega_T \tag{22.7}\]

The initial values of each \(\omega_i\) and of \(\omega_T\) are coupled values. They can be found by solving the implicit, coupled values equations, 22.8 and 22.9. If the design equations include an energy balance on a heat exchange fluid, the exchange fluid energy balance can be added to these coupled values equations they can be solved to additionally find \(T_{ex}\).

\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{i,1} = 0 \tag{22.8}\]

\[ \left(\omega_{i,0},\, \omega_{T,0}\right) \,:\, \omega_{T,1} = 0 \tag{22.9}\]

The Danckwerts boundary conditions,Equations 22.10 and 22.11 can be used to calculate initial values for the concentrations and the temperature. Noting that at stead-state, the mass flow rate of the feed equals the mass flow rate at the reactor inlet leads to Equation 22.12 for the initial value of the volumetric flow rate if the reacting fluid is an ideal gas. If a momentum balance is included among the design equations, the feed pressure can be used as the initial value.

\[ C_{i,0} = \frac{\dot{V}_{feed}C_{i,feed} - D_{ax} \frac{\pi D^2}{4} \omega_{i,0}}{\dot{V}_0} \tag{22.10}\]

\[ 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) } \tag{22.11}\]

\[ \dot{V}_0 = \frac{\dot{m}}{\displaystyle \sum_i \left(M_i C_{i,0}\right)} \tag{22.12}\]

22.3 Axial Dispersion Reactor Model

This section describes the essential components of an axial dispersion reactor model. Upon first reading, it may seem abstract and difficult to follow. If so, it may be helpful to review the examples at the end of the chapter and then re-read this section. Figure 22.2 illustrates the essential components of an axial dispersion reactor model and their relation to numerical equation solvers.

Figure 22.2: Components of an axial dispersion reactor model.

The gray box in Figure 22.2 represents the axial dispersion reactor model. As indicated by the yellow and orange boxes at the bottom of the figure, any constants provided in the assignment narrative, any known universal constants, and any quantities that have been made globally available (see below) can be used anywhere within the reactor model. The three functions included within the model (i. e. the derivatives, coupled value residuals, and reactor functions) must be written by the person solving the design equations. In constast, the two red boxes represent numerical equation solvers. They are not part of the model, but instead are parts of a mathematics software package that is used when solving the design equations. The black arrows indicate that the ATE solver is called by the reactor function and the IVODE solver is called by both the reactor function and the coupled values residuals function. ATE and IVODE solvers are described in Appendices J and K.

The red arrows in Figure 22.2 indicate that the IVODE solver calls the model’s derivatives function and the ATE solver calls the model’s coupled values residuals function. The solvers expect the derivatives function and the residuals function to accept specific arguments and to return the values of specific quantities. Put differently, the only arguments that can be passed to the derivatives function and the residuals function are those specified by the solver mathematics package, and the derivatives and residuals functions must return only the values specified by the solver mathematics package. Appendices J and K identify the arguments that must be passed to these solvers and the values they must return.

As indicated by the dashed orange arrow in Figure 22.2, if any other quantities need to be passed to the derivatives or residuals function, those quantities can be made globally available. Once they are made globally available, they can be used within the derivatives and residuals functions. (There are ways of providing quantities to the derivatives and residuals functions other than by making the quantities globally available, but that is a programming issue and won’t be discussed here.)

To summarize the axial dispersion model can access external given and known constants and globally available quantities. It also interfaces with an external ATE solver and an external IVODE solver. As Figure 22.2 shows, the axial dispersion reactor model has several internal components that now can be described.

Design Equations As with other reactor models, the first component of the axial dispersion reactor model is the design equations. The design equations were listed above, Equations 22.3 through 22.7. If there is pressure drop, a momentum balance must be added to them, as must an energy balance on an exchange fluid that transfers sensible heat at an unknown exchnage fluid temperature. As previously noted, the design equations are a set of CVODEs (see Appendix M). They will be solved using an IVODE solver, so they should be written as derivative expressions (see Appendix K).

Initial Values and Stopping Criterion When solving the CVODE design equations numerically, initial values and a stopping criterion must be provided to the IVODE solver. The initial values for solving the design equations were also listed in the preceding section. They are summarized in Table 22.1 which indicates that the initial values for each \(\omega_i\) and for \(\omega_T\) are coupled values that cannot be calculated directly. If a momentum balance is included in the design equations, the table will have an additional row for the initial value of the pressure. Similarly, if the reacting fluid is a liquid, there will not be a row for the initial value of the volumetric flow rate.

Table 22.1: Initial values and stopping criterion for solving the axial dispersion design equations.

Variable	Initial Value	Stopping Criterion
\(z\)	\(0\)	\(L\)
\(\omega_i\)	\(\omega_{i,0}\) (coupled value)
\(\omega_T\)	\(\omega_{T,0}\) (coupled value)
\(\dot{V}\)	\(\dot{V}_0 = \frac{\dot{m}_{feed}}{\displaystyle \sum_i \left(M_i C_{i,0}\right)}\)
\(C_i\)	\(C_{i,0} = \frac{\dot{V}_{feed}C_{i,feed} - D_{ax} \frac{\pi D^2}{4} \omega_{i,0}}{\dot{V}_0}\)
\(T\)	\(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) }\)

While Table 22.1 indicates the reactor length as the stopping criterion it could also be an outlet concentration or the outlet temperature. The final values of each \(\omega_i\) and of \(\omega_T\) are also known (Equations 22.8 and 22.9) but they cannot be used as the stopping criterion when solving the IVODEs because they are used to solve the coupled value equations.

Derivatives Function In addition to initial values and a stopping criterion, a derivatives function must be provided to the IVODE solver. The derivatives function is essentially the same as for any other reactor type where the design equations are IVODEs. It is used by the IVODE solver. The only arguments passed to it are values of the independent and dependent variables at the start of an integration step, and the only values it returns are those of the derivative of each dependent variable with respect to the independent variable.

In an axial dispersion reactor model, the derivatives function must receive \(z\), each \(C_i\), \(T\), each \(\omega_i\), \(\omega_T\), \(\dot{V}\) (if the system is gas phase), and \(P\) (if there is a momentum balance amoung the design equations). It must return the corresponding values of each \(\frac{dC_i}{dz}\), \(\frac{dT}{dz}\), each \(\frac{d\omega_i}{dz}\), \(\frac{d\omega_T}{dz}\), \(\frac{d\dot{V}}{dz}\) (if \(\dot{V}\) is an argument), and \(\frac{dP}{dz}\) (if \(P\) is an argument). The definition of the derivatives function should include the sequence of calculations that is used to calculate the return values starting from the given and known constants, globally available quantities, and the arguments it receives.

Coupled Values Equations The coupled value equations, 22.8 and 22.9, are reproduced below. They are called coupled value equations because to calculate the \(\omega_{i,0}\) and \(\omega_{T,0}\) values, it is necessary to solve the IVODEs, but the \(\omega_{i,0}\) and \(\omega_{T,0}\) values are needed to solve the IVODEs. The equations are implicit, effectively stating that the desired \(\omega_{i,0}\) and \(\omega_{T,0}\) values are those that, when used to solve the IVODEs, yield values of each \(\omega_{i,1}\) and \(\omega_{T,1}\) that equal zero. They are treated like ATEs and are solved numerically. Hence, they are written as residual expressions.

\[ \left(\omega_{i,0}, \omega_{T,0}\right)\,:\, \omega_{i,L} = 0 = \epsilon_i \]

\[ \left(\omega_{i,0}, \omega_{T,0}\right)\,:\, \omega_{T,L} = 0 = \epsilon_T \]

If the model requires an exchange fluid energy balance involving sensible heat, Equation 22.13, the exchange fluid temperature is also a coupled value. It is coupled because \(T_{ex}\) is needed to solve the IVODEs, but the reactor temperture profile (found by solving the IVODEs) in the first integral below is needed to calculate \(T_{ex}\).

\[ 0 = \int_0^L \pi DU\left( T_{ex} - T \right)dz + \dot m_{ex} \int_{T_{ex,in}}^{T_{ex}} \tilde C_{p,ex}dT = \epsilon_{ex} \tag{22.13}\]

Coupled Values Guess Since an ATE solver is used to solve the coupled values equations, a guess for the coupled values must be provided to it. Specifically guesses are need for each \(\omega_{i,0}\) and for \(\omega_{T,0}\). Looking at Equations 22.6 and 22.7 shows that these variables are the slope of the concentration and temperature profiles at the reactor inlet.

While it is just a guess, if it is too far off the ATE solver may not converge. As was found for steady-state CSTRs, convergence is often most affected by the guess for the temperature. The examples at the end of the chapter illustrate a few strategies for guessing the coupled values.

Coupled Values Residuals Function The other thing that must be provided to the ATE solver when solving the coupled values equations is a residuals function. As is true any time an ATE solver is used, the residuals function must receive only values for the ATE unknowns as arguments, and it must return only the corresponding values of the residuals for the equations being solved. In the case of axial dispersion models, the residuals function will receive values for each \(\omega_i\) and for \(\omega_T\). If the reactor exchanges heat with an external heat exchange fluid, it will also receive the temperature of that fluid as an argument.

Before the coupled values equations can be solved the IVODEs must be solved. Consequently, the residuals function must first calculate each \(C_i\) and \(T\) using Equations 22.10 and 22.11. Then, having all of the initial values, it can call the IVODE solver using those initial guesses along with the stopping criterion and derivatives function as arguments. The IVODE solver will return sets of values of \(z\), each \(C_i\), \(T\), each \(\omega_i\), and \(\omega_T\), spanning the range from \(z=0\) to \(z=L\). If a momentum balance is among the design equations, it also will receive values of \(P\), spanning the range from the reactor inlet to the reactor outlet. The returned values of each \(\omega_i\) and of \(\omega_T\) at the reactor outlet can then be used to evaluate the residuals using Equations 22.8 and 22.9.

Reactor Model Function The reactor model function is the final component of the axial dispersion reactor model. Briefly, it simply needs to solve the coupled values equations using an ATE solver and then use the results to solve the axial dispersion design equations using an IVODE solver. More specifically the reactor model function works as follows.

Guess values for each \(\omega_i\) and for \(\omega_T\).
Get values for each \(\omega_i\) and for \(\omega_T\) by calling an ATE solver, passing the coupled values guesses and the residuals function as arguments.
Calculate the other initial values in Table 22.1.
Get sets of values of \(z\), each \(C_i\), \(T\), each \(\omega_i\), \(\omega_T\), and \(P\) (if there is a momentum balance among the design equations), spanning the range from \(z=0\) to \(z=L\) by calling an IVODE solver passing the initial values, stopping criterion, and derivatives function as arguments.
Return the profiles that the IVODE solver returned.

22.3.1 Multiplicy of Steady States

Just as was seen with steady-state CSTRs, there can be multiple solutions to the steady-state axial dispersion model design equations. Multiple steady states are a consequence of mixing between the reactor inlet and its outlet. They are not observed in steady-state PFRs because there isn’t any axial mixing, but by adding the axial dispersion term, multiple steady states become possible.

22.4 Examples

The first example that follows involves an adiabatic reaction of a liquid feed stream. As such, only the axial dispersion reactor design equations are needed. The second example considers an adiabatic gas phase reaction in a packed bed with a pressure drop. Its solution requires use of the differential ideal gas law and a momentum balance in addition to the axial dispersion reactor design equations.

22.4.1 Liquid Phase Reaction in a Tubular Reactor with Axial Dispersion

22.4.2 Gas Phase Reaction in a Packed Bed Reactor with Axial Dispersion and Pressure Drop

The steady-state feed to a PFR is 75% A and 25% inert gas flowing at 100 cm³ s^-1 , 3 atm and 400 °C. Reaction(1) takes place within the reactor. The reactor diameter is 2.5 cm and its length is 8 m. It is surrounded by a shell containing a molten salt at a constant temperature of 375 °C. The heat transfer coefficient between the reactor and shell is 187 kJ h^-1 m^-2 K^-1. The reaction is catalytic, and the reactor is packed with catalyst particles with a 0.25 cm diameter, a sphericity of 0.7, and a bed porosity of 0.60. The rate expression for reaction (1) is given in equation (2) where \(k_{0,f}\) = 9.0 x 10¹⁷ mol cm³ s^-1 atm^-1, \(E_f\) = 285 kJ mol^-1, \(k_{0,r}\) = 4.09 x 10^-4 mol cm^-3 s^-1 atm^-4, and \(E_r\) = 85 kJ mol^-1. The heat of reaction is constant and equal to 200 kJ mol^-1. The heat capacities of the reagents are also constant: \(\hat{C}_{p,A}\) = 11.7 cal mol^-1 K^-1, \(\hat{C}_{p,Y}\) = 8.3 cal mol^-1 K^-1, \(\hat{C}_{p,Z}\) = 4.3 cal mol^-1 K^-1, and \(\hat{C}_{p,I}\) = 5.8 cal mol^-1 K^-1. The viscosity may be assumed to be constant and equal to 0.027 cP. What are the outlet temperature, pressure and conversion if axial dispersion is included in the model with \(D_{ax}\) equal to ? and \(\lambda_{ax}\) equal to ?? The molecular weights of I, A, Y, and Z are 28, 84, 78, and 2 g/mol, respectively.

\[ A \rightleftarrows Y + 3Z \tag{1} \]

\[ r = k_fP_A - k_rP_YP_Z^3 \tag{2} \]

22.4.2.1 Assignment Summary

Reactor System” Heated steady-State Packed Bed Reactor with Pressure Drop and Axial Dispersion.

Reactor Schematic:

Figure 22.3: Schematic representation of the packed bed reactor.

Quantities of Interest: \(T_1\), \(P_1\), and \(f_A\).

Given and Known Constants: \(y_{A,feed}\) = 0.75, \(y_{I,feed}\) = 0.25, \(\dot{V}_{0}\) = 100 cm³ s^-1, \(\dot{m}_{0}\) = 0.44 g s^-1, \(P_{0}\) = 3 atm, \(T_{feed}\) = 400 °C, \(D\) = 2.5 cm, \(L\) = 8 m, \(T_{ex}\) = 375 °C, \(U\) = 187 kJ h^-1 m^-2 K^-1, \(D_p\) = 0.25 cm, \(\Phi_s\) = 0.7, \(\varepsilon\) = 0.60, \(k_{0,f}\) = 9.0 x 10¹⁷ mol cm³ s^-1 atm^-1, \(𝐸_𝑓\) = 285 kJ mol^-1, \(𝑘_{0,𝑟}\) = 4.09 x 10^-4 mol cm^-3 s^-1 atm^-4, \(𝐸_𝑟\) = 85 kJ mol^-1, \(\Delta H_1\) = 200 kJ mol^-1, \(\hat{C}_{p,A}\) = 11.7 cal mol^-1 K^-1, \(\hat{C}_{p,Y}\) = 8.3 cal mol^-1 K^-1, \(\hat{C}_{p,Z}\) = 4.3 cal mol^-1 K^-1, and \(\hat{C}_{p,I}\) = 5.8 cal mol^-1 K^-1, \(\mu\) = 0.027 cP, \(D_{ax}\) = ?, \(\lambda_{ax}\) = ?, and \(R\) = 8.314 J mol^-1 K^-1 = 82.06 cm³ atm mol^-1 K^1-.

22.4.2.2 Mathematical Formulation of the Calculations

Axial Dispersion Reactor Model

Design Equations

\[ - D_{ax}\frac{d \omega_A}{d z} + \frac{4 C_A}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_A - r \tag{3} \]

\[ - D_{ax}\frac{d \omega_Y}{d z} + \frac{4 C_Y}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_Y + r \tag{4} \]

\[ - D_{ax}\frac{d \omega_Z}{d z} + \frac{4 C_Z}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_Z + 3r \tag{5} \]

\[ - 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 \tag{6} \]

\[ - \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) - r \Delta H_j \tag{7} \]

\[ \frac{d\dot{V}}{dz} = \frac{- \left(M_I \omega_I + M_A \omega_A + M_Y \omega_Y + M_Z \omega_Z \right)}{\left( M_I C_I + M_A C_A + M_Y C_Y + M_Z C_Z\right) } \tag{8} \]

\[ \frac{dC_A}{dz} = \omega_A \tag{9} \]

\[ \frac{dC_Y}{dz} = \omega_Y \tag{10} \]

\[ \frac{dC_Z}{dz} = \omega_Z \tag{11} \]

\[ \frac{dC_I}{dz} = \omega_I \tag{12} \]

\[ \frac{dT}{dz} = \omega_T \tag{13} \]

\[ \frac{{dP}}{{dz}} = - \frac{{1 - \varepsilon }}{{{\varepsilon ^3}}}\frac{{{G^2}}}{{\rho {\Phi _s}{D_p}}}\left[ {\frac{{150\left( {1 - \varepsilon } \right)\mu }}{{{\Phi _s}{D_p}G}} + 1.75} \right] \tag{14} \]

Initial Values and Stopping Criterion

Table 22.2: Initial values and stopping criterion for solving the axial dispersion design equations with \(z=0\) defined as the inlet to the reactor.

Variable	Initial Value	Stopping Criterion
\(z\)	\(0\)	\(L\)
\(\omega_i\)	\(\omega_{i,0}\) (coupled value)
\(\omega_T\)	\(\omega_{T,0}\) (coupled value)
\(\dot{V}\)	\(\dot{V}_0 = \frac{\dot{m}}{M_A C_{A,0} + M_Y C_{Y,0} + M_Z C_{Z,0} + M_I C_{I,0}}\)
\(C_A\)	\(C_{A,0} = \frac{\dot{V}_{0}C_{A,feed} + D_{ax} \frac{\pi D^2}{4} \omega_{A,0}}{\dot{V}_{0}}\)
\(C_Y\)	\(C_{Y,0} = \frac{\dot{V}_{0}C_{Y,feed} + D_{ax} \frac{\pi D^2}{4} \omega_{Y,0}}{\dot{V}_{0}}\)
\(C_Z\)	\(C_{Z,0} = \frac{\dot{V}_{0}C_{Z,feed} + D_{ax} \frac{\pi D^2}{4} \omega_{Z,0}}{\dot{V}_{0}}\)
\(C_I\)	\(C_{I,0} = \frac{\dot{V}_{0}C_{I,feed} + D_{ax} \frac{\pi D^2}{4} \omega_{I,0}}{\dot{V}_{0}}\)
\(T\)	\(T_0 = \frac{\dot{V}_{feed} \left(C_A \hat{C}_{p,A} + C_Y \hat{C}_{p,Y} + C_Z \hat{C}_{p,X} + C_I \hat{C}_{p,I} \right) T_{feed} + \lambda_{ax} \frac{\pi D^2}{4} \omega_{T,0}}{\dot{V}_{feed} \left(C_A \hat{C}_{p,A} + C_Y \hat{C}_{p,Y} + C_Z \hat{C}_{p,X} + C_I \hat{C}_{p,I} \right)}\)
\(P\)	\(P_0\)

Derivatives Function

Arguments: \(z\), \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), \(\omega_T\), \(\dot{V}\), \(C_A\), \(C_Y\), \(C_Z\), \(C_I\), \(T\), and \(P\).

Returns: \(\frac{d\omega_A}{dz}\), \(\frac{d\omega_Y}{dz}\), \(\frac{d\omega_Z}{dz}\), \(\frac{d\omega_I}{dz}\), \(\frac{d\omega_T}{dz}\), \(\frac{d\dot{V}}{dz}\), \(\frac{dC_A}{dz}\), \(\frac{dC_Y}{dz}\), \(\frac{dC_Z}{dz}\), \(\frac{dC_I}{dz}\), \(\frac{dT}{dz}\), and \(\frac{dP}{dz}\).

Algorithm:

\[ \frac{d\dot{V}}{dz} = \frac{- \left(M_I \omega_I + M_A \omega_A + M_Y \omega_Y + M_Z \omega_Z \right)}{\left( M_I C_I + M_A C_A + M_Y C_Y + M_Z C_Z\right) } \tag{8} \]

\[ k_f = k_{0,f} \exp{\left( \frac{-E_f}{RT} \right)} \tag{15} \]

\[ k_r = k_{0,r} \exp{\left( \frac{-E_r}{RT} \right)} \tag{16} \]

\[ P_A = \frac{C_A}{C_A + C_Y + C_Z +C_I}P \tag{17} \]

\[ P_Y = \frac{C_Y}{C_A + C_Y + C_Z +C_I}P \tag{18} \]

\[ P_Z = \frac{C_Z}{C_A + C_Y + C_Z +C_I}P \tag{19} \]

\[ r = k_fP_A - k_rP_YP_Z^3 \tag{2} \]

\[ - D_{ax}\frac{d \omega_A}{d z} + \frac{4 C_A}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_A - r \tag{3} \]

\[ - D_{ax}\frac{d \omega_Y}{d z} + \frac{4 C_Y}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_Y + r \tag{4} \]

\[ - D_{ax}\frac{d \omega_Z}{d z} + \frac{4 C_Z}{\pi D^2}\frac{d \dot{V}}{d z} = - \frac{4 \dot{V}}{\pi D^2}\omega_Z + 3r \tag{5} \]

\[ - 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 \tag{6} \]

\[ - \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) - r \Delta H_j \tag{7} \]

\[ \frac{d\dot{V}}{dz} = \frac{- \left(M_I \omega_I + M_A \omega_A + M_Y \omega_Y + M_Z \omega_Z \right)}{\left( M_I C_I + M_A C_A + M_Y C_Y + M_Z C_Z\right) } \tag{8} \]

\[ \frac{dC_A}{dz} = \omega_A \tag{9} \]

\[ \frac{dC_Y}{dz} = \omega_Y \tag{10} \]

\[ \frac{dC_Z}{dz} = \omega_Z \tag{11} \]

\[ \frac{dC_I}{dz} = \omega_I \tag{12} \]

\[ \frac{dT}{dz} = \omega_T \tag{13} \]

\[ \rho = \sum_i \left(M_i C_i\right) \tag{20} \]

Coupled Values Equations

\[ \left(\omega_{A,0}, \omega_{Y,0}, \omega_{Z,0}, \omega_{I,0}, \omega_{T,0}\right)\,:\, \omega_{A,1} = 0 = \epsilon_A \tag{26} \]

\[ \left(\omega_{A,0}, \omega_{Y,0}, \omega_{Z,0}, \omega_{I,0}, \omega_{T,0}\right)\,:\, \omega_{Y,1} = 0 = \epsilon_Y \tag{27} \]

\[ \left(\omega_{A,0}, \omega_{Y,0}, \omega_{Z,0}, \omega_{I,0}, \omega_{T,0}\right)\,:\, \omega_{Z,1} = 0 = \epsilon_Z \tag{28} \]

\[ \left(\omega_{A,0}, \omega_{Y,0}, \omega_{Z,0}, \omega_{I,0}, \omega_{T,0}\right)\,:\, \omega_{I,1} = 0 = \epsilon_I \tag{29} \]

\[ \left(\omega_{A,0}, \omega_{Y,0}, \omega_{Z,0}, \omega_{I,0}, \omega_{T,0}\right)\,:\, \omega_{T,1} = 0 = \epsilon_T \tag{30} \]

Coupled Values Guess

\[ \omega_{A,0,guess} = -\frac{0.5 C_{A,feed}}{L} \tag{31} \]

\[ \omega_{Y,0,guess} = \frac{0.5 C_{A,feed}}{L} \tag{32} \]

\[ \omega_{Z,0,guess} = 3\frac{0.5 C_{A,feed}}{L} \tag{33} \]

\[ \omega_{I,0,guess} = 0.0 \tag{34} \]

\[ \omega_{T,0,guess} = \frac{10}{L} \tag{35} \]

Coupled Values Residuals Function

Arguments: \(\omega_{A,0}\), \(\omega_{Y,0}\), \(\omega_{Z,0}\), \(\omega_{I,0}\), and \(\omega_{T,0}\).

Returns: \(\epsilon_A\), \(\epsilon_Y\), \(\epsilon_Z\), \(\epsilon_I\), and \(\epsilon_T\).

Algorithm:

Calculate \(C_{A,0}\), \(C_{Y,0}\), \(C_{Z,0}\), \(C_{I,0}\), \(T_0\), and \(\dot{V}_0\), Table 22.2.
Get corresponding sets of values of \(z\), \(C_A\), \(C_Y\), \(C_Z\), \(C_I\), \(T\), \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), \(\omega_T\), \(\dot{V}\), and \(P\) spanning the range from \(z=0\) to \(z=L\) by calling an IVODE solver with the initial values, stopping criterion, and derivatives function as arguments.
Use the returned values of \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), and \(\omega_T\) at \(z=L\) to evaluate the residuals using equations (26) through (30).

Reactor Function

Arguments: none.

Returns: sets of values of \(z\), \(C_A\), \(C_Y\), \(C_Z\), \(C_I\), \(T\), \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), \(\omega_T\), \(\dot{V}\), and \(P\) spanning the range from \(z=0\) to \(z=L\).

Algorithm:

Guess values for the coupled values, e. g. using equations (31) through (35).
Get \(\omega_{A,0}\), \(\omega_{Y,0}\), \(\omega_{Z,0}\), \(\omega_{I,0}\), and \(\omega_{T,0}\) by calling an ATE solver with the guess and the residuals function as arguments.
Calculate the initial values and stopping criterion, Table 22.2.
Get and return sets of values of \(z\), \(C_A\), \(C_Y\), \(C_Z\), \(C_I\), \(T\), \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), \(\omega_T\), \(\dot{V}\), and \(P\) spanning the range from \(z=0\) to \(z=L\) by calling an IVODE solver with the initial values, stopping criterion, and derivatives function as arguments.

Quantities of Interest

Call the reactor function to get sets of values of \(z\), \(C_A\), \(C_Y\), \(C_Z\), \(C_I\), \(T\), \(\omega_A\), \(\omega_Y\), \(\omega_Z\), \(\omega_I\), \(\omega_T\), \(\dot{V}\), and \(P\) spanning the range from \(z=0\) to \(z=L\).
Use the values at \(z=L\) to calculate the quantities of interest.

\[ T_1 = T\big\vert_{z=L} \tag{45} \]

\[ P_1 = P\big\vert_{z=L} \tag{46} \]

\[ f_A = \frac{\dot{V}_0 C_{A,feed} - \dot{V}\big\vert_{z=L} C_A\big\vert_{z=L}}{\dot{V}_0 C_{A,feed}} \tag{47} \]

22.4.2.3 Results, Analysis, and Discussion

The calculations were performed as described above.

22.5 Symbols Used in Chapter 22

Symbol	Meaning
\(t\)	Time.