Compurers and Chrm~ral Enyintwing. Vol. 2. pp.

197-199
0 Pergamon Press Ltd 197W. Printed in Crrat Britam

SHORT NOTE

SOLUTION OF STEADY STATE REACTOR MODELS USING THE

CONTINUATION METHOD AND FEASIBLE STARTING POINT

MORDECHAI SHACHAM
Department of Chemical Engineering, Ben Gurion University of the Negev, Beer-Sheva, Israel

(Received 24 February 1978)

Abstract-The solution of steady state chemical reactor models, using classic numerical methods, (like Newton-
Raphson, Quasi-Newton, etc.), is often very difficult.
It is proposed to select constraints for the independent variables as the first step of the solution, then to use a
combination of the Newton-Raphson and continuation methods, lo obtain the solution, starting from a feasible
initial guess.
Several methods for the solution of steady state reactor models are compared. It is shown that the proposed
method is the most reliable and efficient for locating the solution.

INTRODUCTION first method. It finds only local solutions; it can fail if

Computational experiments with steady state chemical multiple steady states do exist.
reactor models have shown that the solution of such A common limitation of both methods is that they
models, using the classic numerical methods, is very require extremely large computational efforts in function
difficult[5,9, 101. evaluations. This limitation is even greater in a process
Steady state reactor models are represented by a simulation environment, where the reactor model has to
system of nonlinear algebraic equations: be solved several times.
This paper proposes a method for selection of the
f(X) = 0. (1) initial guess in conjunction with the use of the NR
method. For the problems where the NR method fails to
Such a system is frequently obtained by lumping or converge to a physically feasible solution, the use of a
discretization of the reactor models. continuation or imbedding method is investigated.
The best known and most widely-used method for the Several different methods for the solution of a steady
solution of systems of nonlinear equations is the state reactor model are compared by solving three bench
Newton-Raphson (NR) method. mark problems.
The iteration function of the NR method for the solu-
SELECTIONOF THE INITIAL POINT Q
tion of equation (1) is:
Because of the physical nature of the chemical reactor
,&+, =a -X-‘Ra, k =o, 1,2,. (2) model, some or all of the variables in Eq. (I) are subject
to constraints
where k is the iteration number and & is the G deriva-
tive of f at .&, usually termed Jacobian. g(f) > 0. (3)
The NR method converges to the solution f* provided
that & is non-singular for k = 0, 1,2,. . . and that the The existence of such constraints is well known, but
initial approximation .C,is close enough to f* [8]. they are usually not defined until after a solution for (I)
This method is unsatisfactory as a general solution has been obtained. At this stage, the constraints are used
method for reactor models because when starting from to test whether the solution is feasible or not.
an inappropriate initial guess, it tends either to converge If we define the constraints before starting the solu-
to a physically unfeasible solution or not to converge at tion, it will be quite simple to select a feasible initial
all. point. Thus instead of selection of the initial guess,
Two approaches have been used with more success for selection of the constraints has to be considered.
the solution of chemical reactor models. We can distinguish among three different types of
(a) The dynamic model of the reactor is integrated constraints according to their origin:
until steady state is obtained. Using this approach, the (I) Those originating from physical and chemical
dynamic operation of the reactor is simulated. Assuming principles.
that the model is sound and that steady state solution (2) Those originating from the limitations of the
does exist, this method will most often obtain the solu- mathematical model.
tion. Ham[5], Smith et al.[ lo] and Silverberg[9] found (3) Extremum or limiting values of the variables for
this method to be the most reliable method for chemical the specified conditions.
reactor simulation. Examples of constraints of the first type are well
(b) The steady state reactor problem is converted into known. The mole fraction of a component cannot be
a minimization problem, which is then solved by a search smaller than zero or larger than unity, negative vapor-
method. This method seems to be less reliable than the liquid equilibrium ratio is meaningless, etc.
197
198 MORDECHAf SHACHAM

An example of constraints of the second type is that dynamic operation of the chemical reactor. The
we do not take into consideration abnormal operation of conclusion from this analogy is that if the initial ap-
a water-cooled reactor, where the cooling water starts to proximation 5 is located inside the feasible region, the
boil. Thus, the outlet temperature of the cooling water continuous solution curve of Eq. (5) cannot leave this
should be below its boiling point temperature. region and the steady state solution for 8” = 0 is
The third type of constraints is best represented by an obtained inside the region.
example. Consider a well-mixed reactor where the reac- While it is diicult to prove this analogy formally, the
tion computational experiments show that it is sound.
According to the proposed algorithm, an initial guess
A+B=CtD .& whi_ch satisfies the constraints (Eq. 2) is specified.
Then Y0 is calculated and the series of problems defined
occurs. Provided that the feed to the reactor does not in Eq. (5) are solved; the first time using At9= 1 (which is
contain any component C, the mole fraction of C in the simply the NR method). If no feasible solution is found
outlet stream cannot exceed the value of 0.5. Using these this way or a singular point has been attained, then A@is
three types of constraints, fairly close boundaries for the divided by ten and the solution procedure starting from
feasible region can be defined. %, is repeated.
For cases where even for the smallest value of A0 an
THE CO~~UATION METHOD unfeasible solution is obtained, a new starting point
If the NR method fails to converge to feasible solution which satisfies the constraints is generated. (This last
even from a feasible starting point, the use of the step of the algo~thm is added only for the sake of
continuation or imbedding methods should be completeness; it was not required in any of the test
considered. The use of these methods for the solution of problems.)
a system of non-linear equations is quite recent. Growing
interest in these methods has been shown since l%O. COZEN OFTHE SOLUTIONMETHODS
Using the continuation method, the system of alge- Five solution methods have been compared: (1) NR
braic Eqs. (1) is converted into a system of differential method, (2) the continuation method, (3) numerical in-
equations. tegration of the dynamic model, (4) Davidon’s method
and (5) the complex method.
The first two have already been discussed and a short
description of the others follows.
If the chemical reactor problem can be considered as
Equ$ion (4) is “integrated” in a special way. For each the steady state solution of a dynamic model, then the
step (&/dt)’ is spe&ed (i-is the number of the in- dynamic model [which is similar to Eq. (4)] can be
tegration step), the specified value being: numerically integrated until steady state is achieved. The
CSMP simulation program was used for the integration
of the equation of the dynamic model.
Equation (1) can be converted into the following
where 8’ = 1; 8’ = 8’~’ -A& and 8” = 0, minimization probIem:
Equation (4) is solved for f’ using the NR method. The
sequence of the problems which have to be solved is Minimize S’ = Cfi(Q’ (8)
shown in Eq. (5).
subject to the constraints given in Eq. (3).
c?(f,e’)=f(f’f-eiFO=o i=l,2 ,..., n. 6% The min~um of S* is clearfy the solution of Eq. (1).
This minimization problem can be solved using non-
For 8” =0 the roots of Eqs. (5) and (1) are identical. linear opt~~tion techniques. The FMFP subroutine of
This system is solved for increasing values of i. The first IBM[7] and the OPTIM program of EvansD] were used
approximation for each step: & is obtained from the for the soIution of the migration problem. The FMFP
solution of the earlier step. subroutine uses the unconstrained minimization method
of Davidon which was described by Fletcher &
ti = f*.i- I PowellI41. OPTIM uses the Complex method of Boxfl]
(6)
which is a constrained search technique.
or, if i 2 2 by linear extrapolation: Three problems have been solved using the above
methods. Example No. 1 is an isothermal backmix reac-
&o’= 2z*.‘- I _ j*.i-z tor of the Williams-Otto process, taken from p. 408 of
(7)
[6]. In this problem there are six variables, all of which
It can be shown (181,p. 336) that if AB is chosen small represent reactor outlet Aow rate. The constraints for the
enough, the continuation method, as defined by Eqs. (5) variables as well as the solution of all the examples are
and (6), always converges regardless of the initial ap- shown in Table 1,
proximation. Thus, one of the restrictions of the NR Example No. 2 is a backmix, jacketed reactor in which
method is overcome by the continuation method. an exothermal reaction takes place. The lirst four vari-
However singular points may still present problems. ables in this example represent component liow rates, the
From the above, it is clear that the continuation last two variables are the temperature inside the reactor
method has better convergence properties than the NR and the cooling water temperature, respectively. This
method. It must be shown now that the continuation example is taken from [lo].
method has to converge to a feasible solution. Example No. 3 deals with partial oxidation of
As can be seen from Eq. (4), the continuation method methane. It is required to find the O&H4 reactant ratio
may be considered approx~ately as a s~ulation of the that will produce an adiabatic equ~ib~um tem~rature of
 Solution of steady state reactor models 199

Table I. Constraints and solution values of the variables for the problems

Example 1 Example 2 Example 3

Variable Lower Upper Lower Upper Lower Upper
No. limit Solution limit limit Solution limit limit Solution limit

I 0 16.22 19.9 0 2.356 3.6 0 0.3229 0.33

2 0 44.17 51 0.4 0.992 2.2 0 0.092 0.5
3 0 2.97 70 0 0.609 1.8 0 0.046 0.5
4 40 46.8 75 0.48 0.486 0.49 0 0.618 0.67
5 0 1.27 35 7s 106.9 220 0 0.037 0.67
h 4 6.68 40 75 170.8 318 0.5 0.577 2
? 1.5 2.9778 3

Table 2. Results for Examples 1and 2

Method of solution NR CSMP FMFP OPTIM

-
Example No. I No. of iterations* 5 I50 540
CP time (s)* 0.07 2.0 3.2 3.5
Example No. 2 No. of iterations* 5 500
CP time (s)* 0.06 2.8 4.2
-
*Average values for several runs starting from different initial points. The
central processor (CP) time is based on a CDC Cyber 73 computer.

2200°F. (This example is taken from p. 321 of [2].) There best is the continuation method. The NR method, FMFP
are seven variables; the first five are mole fractions of and OPTIM converge to the correct solution very rarely
different components and the last two represent mole under these circumstances.
ratios of some of the components. This example is It may be concluded that the selection of the con-
different from the first two, since this problem cannot be straints for the variables is an essential step in the
represented by a dynamic model. solution of steady state chemical reactor problems. Pro-
vided that the initial estimate is located inside the feas-
RESULTS AND DlSCUSSlON ible region and that the problem does have a solution,
The three examples have been solved severai times the use of combined NR and continuation methods is the
using different but feasible starting points. most reliable and efficient manner for locating the solu-
The results for Examples 1 and 2 are shown in Table 2. tion
Since the NR method converged to the correct solution
in ali of these cases, there was no need to use the REFXRENCES
continuation method. 1. M. J. Box, D. Davies & W. H. Swann, Non-Linear Opti-
The other three methods also converged to the correct misation Techniques, I.C.1. Monograph No. 5. Oliver &
solution, but they used up 20-70 times more CP time Boyd, Edinburgh (1%9).
than the NR method. 2. B. Camahan, H. A. Luther & 1. 0. Wilkes, Applied Numeti-
Example No. 3 proved to be the most difficult of the cal Methods. John Wiley,New York (1%9).
three and the dynamic model method could not be used 3. L. B. Evans, Private communications(1975).
4. R. Fletcher & M. J. D. Powell, Comput. Jl 6(2), 163-168
for this problem. FMFP and OPTIM did not converge in
(l%3).
1000 iterations regardless of the starting point used. The 5. P. G. Ham, The transient analysis of integrated chemical
continuation method converged from any feasible start- processes, Ph.D. Dissertation, Univ. Pennsylvania,
ing point in about 30 iterations (0.32CP s). The NR Phil~elphia (t971).
method converged only from a starting point which was 6. E. J. Henley &IE. M. Rosen, Material and Energy Balance
close to the mid~int between the lower and upper limits C~~putat~n. John Wiley, New York (1969).
(6 iterations, 0.08 CP s). Starting the NR method from 7. I.B.M. System/360 Scientific Subroutine Package, Form No.
initial points close to the lower or upper limit lead to GH’20-0205-4, 5th Ed. (1970).
convergence to a singular point or unfeasible solution. 8. J. M. Ortega & W. C. Rheinboldt, iterutioe Solution of
Nonlinear Equations in Several Variables. Academic Press,
Some experiments were performed to start with un-
New York (1970).
feasible initial estimations. In general, the different 9. A. Silverberg, Development of a general method for the
methods converged to the correct solution in some cases design of a chemical reactor using a computer, MSc. Thesis
and diverged or converged to unfeasible solutions in fin Hebrew),Technion, Haifa (1975).
others. If we compare the different methods on the basis 10. C. L. Smith, R. W. Pike & P. W. Murrill, Formulation and
of convergence from unfeasible points, then the dynamic ~ptimizution of ~uthematjcai yodels. Internation~ Text-
model approach proves to be the best and the second book Co., Scranton (IWO).