Appendix I — Solving Algebraic/Transcendental Equations

This appendix examines the numerical solution of sets of algebraic/transcendental equations (ATEs). There are numerous software packages that include an ATE solver (function that solves a set of ATEs numerically). While the details differ from one software package to another, the vast majority of ATE solvers require the same input, use it the same way, and return equivalent results. In order to follow Reaction Engineering Basics example assignments wherein sets of ATEs ae solved, and in order to implement them using an ATE solver of one’s choosing, a general understanding of how ATE solvers work and the input that must be provided to them is needed. The intent of this appendix is to provide that necessary understanding.

I.1 Identifying ATEs

In Reaction Engineering Basics, almost any equation that does not contain an ordinary derivative or a partial derivative is likely an ATE. Slightly more specifically, a set of ATEs is a group of 1 or more mathematical equations that may involve or contain math operations (addition, subtraction, multiplication, and division), quantities raised to powers, and transcendental functions. Exponential functions are the most common transcendental functions appearing in the ATEs found in Reaction Engineering Basics. They arise any time an equation includes a rate coefficient that displays Arrhenius temperature dependence, Equation 4.8.

I.2 Defining Residuals for ATEs and a Residuals Function

When solving a set of ATEs numerically, it is useful to define a residual that corresponds to each ATE. To do so, rearrange the ATE so that there is a zero on one side of the equation. The residual is equal to the non-zero side of the resulting equation. A residual can be defined for each ATE in the set to be solved. Note that in general each residual is a function of all of the unknowns. A guess for a solution of the ATEs can be tested by using the guess to calculate the residuals. If all of the residuals equal zero, the guess is a solution of the ATEs.

A computer function, here called a residuals function, can be written and used to test guesses in that way. The residuals function receives as an argument, a guess for the solution. Within the residuals function, the guess is used to calculate the residual corresponding to each of the ATEs, and the resulting values are returned. Often the guesses are passed to the residuals function as a vector and the corresponding residuals are returned as a vector.

I.3 ATE Solvers

Typically, an ATE solver must be provided with two inputs. The first is an initial guess for the unknowns. The second is a residuals function, as defined above. The manner in which these inputs are provided to the solver depends upon the specific solver being used.

Effectively, an ATE solver finds the solution of a set of ATEs by trial and error.

• It tests the initial guess by calling the residuals function and determining whether all of the residuals equal zero.
• If the initial guess is not a solution of the ATEs, the ATE solver generates an improved guess and tests that.
  • The improved guess will not be an exact solution (unless the ATEs are linear equations), but usually it will be closer to the solution than the previous guess.
• The ATE solver keeps generating improved guesses and testing them until no further improvement is possible.
• It then returns the final, most-improved guess as the best solution of the ATEs.

So, the numerical solution of a set of ATEs is an iterative process, and ideally, the residuals should get closer and closer to zero with each iteration. This is called convergence to a solution. Typically, the solution returned by the solver will not be exact, but it will be “very close” to the exact solution. Put differently, when a converged solution is found, the difference between the solution returned by the solver and the exact solution is negligible.

I.3.1 Solver Issues

Sometimes the solver is unable to converge to the point where the residuals are “very close” to zero. In this situation, the solver eventually has to quit without finding a converged solution. Typically it will print an error message or return a flag variable, and that variable will indicate whether or not a converged solution was obtained.

When the ATE solver is not able to converge to a solution, it is often because the guess provided to the solver wasn’t close enough to the solution or because of a singularity in one or more of the residuals. For example, fractions that contain one or more unknowns in the denominator may appear in the residuals. If a guess for the solution causes such a denominator to equal zero, the residual becomes infinite.

Potential singularities can be removed by re-defining the residuals. To do so, first rearrange each ATE so there is a zero on one side of the equation. Then multiply through the equation by each denominator that contains an unknown. Finally, use the resulting non-zero side of the equation as the residual. If possible, a better initial guess should be defined, and the solver should be called again.

Additionally, sets of nonlinear equations can have more than one solution. For example, a quadratic polynomial has two solutions. If a set of ATEs has multiple solutions, most ATE solvers will only find one of them. To find other solutions, the solver should be re-started using a different initial guess.

I.4 Formulating the Solution of a Set of ATEs Mathematically

In general, a set of \(N\) ATEs can contain any number of known, constant quantities. The set of ATEs can be solved to find the values of \(N\) unknown quantities (the unknowns). To do so numerically, it will be necessary to evaluate the residuals given a guess for the solution. Starting with a set of \(N\) residuals corresponding to the \(N\) ATEs the mathematical formulation of the solution proceeds as follows:

1. Identify all known constants and their values.
2. Identify the \(N\) unknowns to be found by solving the \(N\) ATEs.
3. For each quantity appearing in the ATEs other than known constants and the \(N\) unknowns, write an ancillary equation that can be used to calculate it using the known constants and the \(N\) unknowns.
4. Describe how to use the ancillary equations to evaluate the residuals.

I.5 Implementing the Solution of a Set of ATEs Numerically

Numerical implementation of the solution involves writing computer code that performs the calculations as described in the mathematical formulation. That code will work together with an ATE solver. The specific details for writing it will depend upon the programming language and the specific ATE solver being used. Here the implementation is presented in general terms that provide sufficient detail so that the solution can be implemented using mathematics software of the reader’s choosing. Reaction Engineering Basics assumes that the code will be broken into computer functions that each perform a specific task.

reactor model function

sets initial guesses for the unknowns; calls the ATE solver, passing the initial guesses and the name of the residuals function; checks that the ATE solver converged; and returns a solution. That is, it returns one set values of the unknowns that solve the ATEs.

residuals function

receives guesses for the unknowns, uses them to evaluate the residuals, and returns the residuals.

top-level code (and/or a master calculations function)

establishes a mechanism for making all known constants available to all functions, calls the reactor model function to get the solution of the ATEs, and performs any additional calculations that are necessary.

I.6 Example: Solving 3 ATEs to Find 3 Unknowns

Equations (1), (2), and (3) are representative of a set of ATEs that an engineer might generate and need to solve when completing a reaction engineering assignment. For this example, suppose that the engineer needs to find the values of \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\) by solving the ATEs. For this example we will also assume that the engineer knows the rate expression for \(r\) that is used in the solution that follows and the constant values of the following quantities: \(\dot n_{A,in}\) = 500 mol h^-1, \(V\) = 500 L, \(\dot n_{Z,in}\) = 0 mol h^-1, \(\breve{C}_p\) = 1170 cal L^-1 K^-1, \(\dot{V}\) = 250 L h^-1, \(T_{in}\) = 423 K, \(\Delta H\) = 18,200 cal mol^-1, \(k_0\) = 1.14 x 10⁹ L mol^-1 h^-1, and \(E\) = 16,200 cal mol^-1.

\[ \dot n_{A,in} = \dot n_A + Vr \tag{1} \]

\[ \dot n_{Z,in} = \dot n_Z - Vr \tag{2} \]

\[ \breve{C}_p \dot{V}\left( T - T_{in} \right) = - Vr\Delta H \tag{3} \]

---

I.6.1 Assignment Summary

Click Here to See What an Expert Might be Thinking at this Point

The narrative for this example assigned appropriate variable symbols to every quantity it provided. I can see that I’m going to need the ideal gas constant, so I’ll add it to the given constant quantities. The narrative doesn’t say so, but the equations are the reactor design equations for a steady-state CSTR.

Given and Known Constants: \(\dot n_{A,in}\) = 500 mol h^-1, \(V\) = 500 L, \(\dot n_{Z,in}\) = 0 mol h^-1, \(\breve{C}_p\) = 1170 cal L^-1 K^-1, \(\dot{V}\) = 250 L h^-1, \(T_{in}\) = 423 K, \(\Delta H\) = 18,200 cal mol^-1, \(k_0\) = 1.14 x 10⁹ L mol^-1 h^-1, \(E\) = 16,200 cal mol^-1, and \(R\) = 1.987 cal mol^-1 K^-1.

Reactor System: Steady-state CSTR

Quantities of Interest: \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\)

I.6.2 Mathematical Formulation of the Solution

Click Here to See What an Expert Might be Thinking at this Point

Normally I would need to generate the reactor design equations at this point, but here the design equations are given. I will, however, rewrite them as residuals since I know I’ll need them in that form when I solve them numerically

Reactor Design EQuations

Rewriting equations (1), (2), and (3), as residuals expressions yields equations (4), (5), and (6). They will be solved to find the values of \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\).

\[ 0 = \dot n_{A,in} - \dot n_A - Vr = \epsilon_1 \tag{4} \]

\[ 0 = \dot n_{Z,in} - \dot n_Z + Vr = \epsilon_2 \tag{5} \]

\[ 0 = \breve{C}_p \dot{V}\left( T - T_{in} \right) + Vr\Delta H = \epsilon_3 \tag{6} \]

Click Here to See What an Expert Might be Thinking at this Point

In order to solve the ATEs numerically I’ll need to do two things: calculate the values of the residuals for each guess and provide an initial guess.

Starting with evaluating the residuals, I’ll need to evaluate them each time a new guess for the solution is generated. At that point, guesses for each of the unknowns will be available along with the given and known constants. I will need to calculate any other quantities that appear in the residuals expressions.

In this example, the only other quantity that appears in any of the ATEs is \(r\), so I will write an expression for \(r\) in terms of the known constants and the unknowns. (As mentioned in the problem statement, it is assumed here that the engineer knows the rate expression used in equation (7) below.)

Writing the ancillary equation for \(r\) introduces \(C_A\), and \(C_A\) is neither a known constant nor one of the unknowns, so an expression is needed for \(C_A\), too. In this case, the defining equation for concentration in a closed system can be used.

Ancillary Equations for Evaluating the Residuals

\[ r = k_0 \exp{\left( \frac{-E}{RT} \right)} C_A^2 \tag{7} \]

\[ C_A = \frac{\dot{n}_A}{\dot{V}} \tag{8} \]

Click Here to See What an Expert Might be Thinking at this Point

When I solve the ATEs numerically, I’ll also need to provide an initial guess for the solution. I’ll guess that half of the A reacts. That will result in an amount of Z equal to half of the original amount of A. I don’t have any idea how large the temperature change will be, but having some knowledge of the CSTR design equations, I can see that this reactor is adiabatic and the reaction is endothermic. That means that the temperature will decrease. I’ll gues the temperature change is small, 10 K.

Notes

1. Reactor models can be very sensitive to temperature due to the exponential term. For this reason it is prudent when guessing a temperature to guess a value that is close to a known temperature in the system, such as \(T_{in}\) in this example. If the guess is too large, the exponential term may become excessively large and prevent the solver from converging.
2. In this situation (an endothermic reaction in an adiabatic CSTR) the engineer would not expect multiple solutions to be possible. If the engineer knew that multiple solutions were possible, the solution would be repeated using different initial guesses to find additional solutions.

Calculations Summary

1. Substitute the known constants into all equations.
2. When it is necessary to evaluate the residuals
   1. Calculate \(C_A\) using equation (8),
   2. Calculate \(r\) using equation (7), and
   3. Calculate \(\epsilon_1\), \(\epsilon_2\), and \(\epsilon_3\) using equations (4), (5), and (6).

I.6.3 Numerical Implementation of the Solution

1. Make the given and known constants available for use in all functions.
2. Write a residuals function that
   1. receives guess values for \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\),
   2. evaluates the residuals as described in step 2 of the calculations summary, and
   3. returns the values of the residuals.
3. Write a reactor model function that
   1. defines an initial guess for the unknowns, \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\),
   2. gets a set of values of the unknowns that solve the ATEs by calling an ATE solver and passing the following information to it
      1. the initial guess and
      2. the name of the residuals function from step 2 above, d, checks that the solver converged, and
   3. returns the values returned by the ATE solver.
4. Perform the analysis by
   1. getting \(\dot{n}_A\), \(\dot{n}_Z\), and \(T\), by calling the reactor model function from step 3 above,
   2. displaying and saving the results, as desired.

I.6.4 Results and Discussion

The calculations were performed as described above and it was found that \(\dot{n}_A\) = 157 mol/h, \(\dot{n}_Z\) = 343 mol/h, and \(T\) = 128 °C.

Note that the solution is not exact, but the error is negligibly small. That is, if the solution is used to calculate the residuals using equations (4) through (6), the values are not exactly zero, but instead (with the software used here), \(\epsilon_1\) = -4.32 x 10^-10, \(\epsilon_2\) = 4.32 x 10^-10, \(\epsilon_3\) = 7.87 x 10^-6.

SCoRE Connection

Videos showing how to complete this assignment using either Matlab or Python, along with the Matlab and Python code, are available in SCoRE

I.7 Symbols Used in Appendix I

Symbol	Meaning
\(k_0\)	Pre-exponential factor in the Arrhenius expression.
\(\dot n_i\)	Molar flow rate of reagent \(i\); an additional subscripted \(in\) denotes the inlet molar flow rate.
\(r\)	Rate of reaction.
\(\breve{C}_p\)	Volumetric heat capacity of a fluid.
\(C_i\)	Molar concentration of reagent \(i\).
\(E\)	Activation energy in the Arrhenius expression.
\(N\)	Number of ATEs being solved and number of unknowns being found.
\(R\)	Ideal gas constant.
\(T\)	Temperature; an additional subscripted \(in\) denotes the inlet temperature.
\(V\)	Volume of reacting fluid.
\(\dot V\)	Volumetric flow rate.
\(\epsilon_i\)	Residual for the \(i^{th}\) ATE.
\(\Delta H\)	Heat of reaction.