Appendix M — Solving Coupled-Value ODEs

This appendix examines the numerical solution of sets of IVODEs wherein both the initial and the final value of one or more dependent variables are unknown, but are coupled via an ATE. In Reaction Engineering Basics equations of this type are referred to as coupled value ordinary differential equations (CVODEs), and they are only encountered when thermally back-mixed and recycle PFRs. This appendix shows that CVODEs can be solved numerically using an ATE solver (see Appendix J) in conjunction with an IVODE solver (see Appendix K). No other solvers are needed.

M.1 Identifying CVODEs

The sets of CVODEs that are solved in Reaction Engineering Basics are linear, first-order ODEs, as was the case for IVODEs in Appendix K. CVODEs are encountered when modeling thermally back-mixed PFRs and Recycle PFRs. The distinguishing feature of CVODEs is that for at least one of the dependent variables, neither the initial nor the final value is known or can be calculated directly, but they are coupled through an ATE.

In a thermally back-mixed PFR, the initial and final values of the temperature are not known and cannot be calculated directly, but they are coupled through an energy balance on the heat exchanger in the system (see Chapter 16). In a recycle PFR the initial and final values of the molar flow rates and the temperature are coupled through mole and energy balances on the mixing point in the sustem (see Chapter 17).

Both the initial values and the final values appear in the ATEs that couple them. No special preparation is required when solving CVODEs. The key to solving them is to choose the IVODE initial values as the unknowns in the ATEs. By doing so, values for the IVODE initial values will be provided to the residuals function as arguments. Using those values the IVODE design equations can be solved within the residuals function and the IVODE final values can be extracted from the results. Then the ATE residuals can be evaluated.

M.2 Mathematical Formulation of the Solution of a Set of CVODEs

When solving CVODEs, the solution of the ATEs is formulated just as described in Appendix J and the solution of the IVODEs is formulated as described in Appendix K. The first important aspect of the formulation is that the IVODE initial values must be passed as arguments to the function that solves the IVODEs. The second important aspect is that by choosing the IVODE initial values as the ATE unknowns, values for the ATE unknowns will be passed to the residuals function as arguments. The first step in the residuals function algorithm then is to solve the IVODEs to find their final values. This facilitates evaluation of the residuals. Once the IVODE initial values have be found by solving the ATEs, they can be used to solve the IVODEs.

M.3 Example

The following example illustrates the mathematical formulation of the solution of a set of CVODEs as described above.

M.3.1 Formulation of the Numerical Solution of Recycle PFR Design Equations

A reaction engineer needs to solve the recycle PFR design equations shown in equations (1) through (3) using the initial values and stopping criterion in Table M.1. The initial values of the dependent variables are related to their final values through the mixing point balances shown in equations (4) - (6). The rate can be calculated using equation (7). The values of \(D\), \(\Delta H\), \(R_R\), \(\breve{C}_p\), \(k_0\), \(E\), \(R\), \(L\), \(\dot{V}_{feed}\), \(\dot{n}_{A,feed}\), and \(T_{feed}\) are known constants.

\[ \frac{d\dot{n}_A}{dz} = - \frac{\pi D^2}{4}r \tag{1} \]

\[ \frac{d\dot{n}_Z}{dz} = \frac{\pi D^2}{4}r \tag{2} \]

\[ \frac{dT}{dz} = - \frac{\pi D^2 r \Delta H}{4 \left(1+R_R\right) \dot{V}_{feed} \breve{C}_p} \tag{3} \]

Table M.1: Initial values and stopping criterion for solving equations, (1) - (3).
Variable Initial Value Stopping Criterion
\(z\) \(0\) \(L\)
\(\dot{n}_A\) \(\dot{n}_{A,0}\)
\(\dot{n}_Z\) \(\dot{n}_{Z,0}\)
\(T\) \(T_0\)

\[ 0 = \dot{n}_{A,0} - \dot{n}_{A,feed} - \frac{R_R}{1 + R_R} \dot{n}_{A,1} = \epsilon_1 \tag{4} \]

\[ 0 = \dot{n}_{Z,0} - \frac{R_R}{1 + R_R} \dot{n}_{Z,1} = \epsilon_2 \tag{5} \]

\[ T_0 - T_{feed} + R_R \left( T_0 - T_1 \right) = \epsilon_3 \tag{6} \]

\[ r = k_0 \exp{ \left( \frac{-E}{RT} \right)} \frac{\dot{n}_A}{\dot{V}_{feed}\left(1+R_R\right)} \tag{7} \]

The IVODE initial values, \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\), and the IVODE final values, \(\dot{n}_{A,1}\), \(\dot{n}_{Z,1}\), and \(T_1\), are unknown, but they are coupled through the mixing point ATEs, equations (4) through (6), making equations (1) through (3) CVODEs. Formulate their solution to find the IVODE initial and final values

PFR Design Equations

Required input: \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\)

Derivatives Function

\(\qquad\) Arguments: \(z\), \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\).

\(\qquad\) Return values: \(\frac{d\dot{n}_A}{dz}\), \(\frac{d\dot{n}_Z}{dz}\), and \(\frac{dT}{dz}\).

\(\qquad\) Algorithm:

  1. Calculate the unknowns in the ODEs

\[ r = k_0 \exp{ \left( \frac{-E}{RT} \right)} \frac{\dot{n}_A}{\dot{V}_{feed}\left(1+R_R\right)} \tag{7} \]

  1. Evaluate and return the derivatives using equations (1) - (3).

Function that solves the ODE design equations

\(\qquad\) Argument: \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\).

\(\qquad\) Algorithm:

  1. Call an IVODE solver
    1. Passing the initial values, stopping criterion, and the derivatives function as arguments.
    2. Receiving \(\underline{z}\), \(\underline{\dot{n}}_A\), \(\underline{\dot{n}}_Z\), and \(\underline{T}\) spanning the range from \(z=0\) to \(z=L\).
  2. Return the results from the IVODE solver.

Mixing Point Balances

Residuals Function

\(\qquad\) Arguments: \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\).

\(\qquad\) Return values: \(\epsilon_1\), \(\epsilon_2\), and \(\epsilon_3\).

\(\qquad\) Algorithm:

  1. Call the function that solves the PFR design equations, above
    1. passing \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\) as the argument
    2. receiving \(\underline{z}\), \(\underline{\dot{n}}_A\), \(\underline{\dot{n}}_Z\), and \(\underline{T}\)
  2. Extract the IVODE final values.

\[ \dot{n}_{A,1} = \dot{n}_A \big \vert_{z=L} \tag{8} \]

\[ \dot{n}_{Z,1} = \dot{n}_Z \big \vert_{z=L} \tag{9} \]

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

  1. evaluate and return \(\epsilon_1\), \(\epsilon_2\), and \(\epsilon_3\) using equations (4) - (6).

Calculating the Quantities of Interest

  1. Guess values for \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\).
  2. Call an ATE solver passing the guess and the residuals function, above, as arguments and receiving \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\).
  3. Call the function, above, that solves the PFR design equations passing \(\dot{n}_{A,0}\), \(\dot{n}_{Z,0}\), and \(T_0\) as the argument and receiving \(\underline{z}\), \(\underline{\dot{n}}_A\), \(\underline{\dot{n}}_Z\), and \(\underline{T}\).
  4. Extract \(\dot{n}_{A,1}\), \(\dot{n}_{Z,1}\), and \(T_1\) using equations (8) - (10).

M.4 Symbols Used in this Appendix

Symbol Meaning
\(k\) Rate coefficient.
\(z\) Axial distance from the reactor inlet.
\(C_A\) Concentration of reagent A; an additional subscripted \(in\) denotes the concentration at the reactor inlet.
\(L\) Length of the reactor.