Singh, C. P. P., and D. N. Saraf, “Simulation of High-Temperature Wen, C. Y.

, “Noncatalytic Heterogeneous Solid Fluid Reaction Models,”

Water-Gas Shift Reactors,”lnd. Eng. Chem., Process Des. Dev., 16, No. Znd. Eng. Chem., 60, No. 9, 34 (1968).
3,313 (1977). Wen, C. Y., and E. S. Lee, Eds., Coal Conversion Technology, Addison
Smoot, L. D., and P. J. Smith, “Mixing and Gasification of Coal in Entrained Wesley Publishing Go. (1979).
Flow Systems,” ERDA Rpt., FE-2666-f, 2, (1979). Wen, C. Y., and L. T. Fan, Models for Flow Systems and Chemical Re-
Sprouse, K. M., “Modeling Pulverized Coal Conversion in Entrained actors, Marcel Dekker, New York (1975).
Flows,” AZChE j . , 26,964 (1980). Wen, C. Y., and T. 2. Chuang, “Entrainment Coal Gasification Modeling,”
Szekely, J., J. W. Evans, and H. Y. Sohn, Gas-Solid Reactions, Academic Ind. Eng. Chem. Process Des. Dev., 18,684 (1979).
Press (1976). Zahradnik, R. L., and R. J. Grace, “Chemistry & Physicsof Entrained Coal
Thring, N. W., and R. H. Essenhigh, Chemistry of Coal Utilization-Sup- Gasification,” Adv. Chem. Ser., 131, 126 (1974).
plementary Volume, H. H. Lowry, Ed., Chap. 17, John Wiley & Sons
(1963).
Ubhayaker, S. K., D. B. Sticker, and R. E. Gannon, “Modeling of Entrained Manuscript received February 17,1982; revision received January 7 , and accepted
Bed Pulverized Coal Gasifiers,” Fuel, 56,281 (1977). February 28,1983.

Decomposition of Systems of Nonlinear

Algebraic Equations
A new method for decomposing irreducible subsets, in the solution of systems MORDECHAI SHACHAM
of nonlinear algebraic equations, is presented. This method consists of two steps:
(1) elimination of the nonlinearity from some of the equations by replacing non- Department of Chemical Engineering
linear expressions by new variables; and (2) formulation of a problem of smaller Ben Gurion University of the Negev
dimension by tearing the linear subset of equations. It is shown that these modi- Beer Sheva, Israel
fications do not change considerably the convergence rate of the Newton-Raphson
and Broyden’s methods while reducing the problem’s dimension. Computer time
reduction up to 80% is demonstrated in the examples solved. An algorithm for
elimination of nonlinear expressions, which uses Boolean matrices instead of
formula manipulation, is also presented.

SCOPE
Steady-state simulation or design of chemical processes gives partitioning.
rise to the need of solving large sets of nonlinear algebraic Irreducible subsets can further be decomposed by “tearing.”
equations. In the equation-oriented flowsheeting programs This operation reduces essentially the dimension of the irre-
which are currently being developed (Shacham et al., 1981), ducible subset, by selecting some of the variables as “tear
these may be problems containing tens of thousand of linear and variables” and expressing the additional variables in the subset
nonlinear equations. The simultaneous solution of such a large as explicit functions of the tear variables. There are several
number of equations is very time-consuming and may be even difficulties in applying tearing techniques to nonlinear systems,
impractical. Fortunately, large systems of equations tend to be the major ones being: (1)tearing requires formula manipulation
sparse so that each equation usually involves a small number which is very expensive in terms of computer time; and (2)the
of variables. This sparsity makes it possible to decompose the smaller system generated by tearing tends to be more nonlinear
large system into a set of smaller problems which can then be and more difficult to solve (Mah, 1972).
solved sequentially. This process (usually denoted partitioning) In this paper a new tearing method for systems of nonlinear
decomposes the system into “irreducible subsets” which cannot equations is presented. This method does not require formula
be further partitioned. It has been proved that it is advantageous manipulation, and it either improves or does not change the
to partition a system to irreducible subsets (Hernandez and convergence rate when the Newton-Raphson or Broyden’s
Sargent, 1979), and there are several algorithms available for methods are used for solution.

CONCLUSIONS AND SIGNIFICANCE

The proposed tearing method includes the following steps: using partial Gauss elimination in this subset to express all the
1. Some of the nonlinear equations are linearized by re- variables as explicit linear functions of the tear variables. The
placing nonlinear expressions by new variables, and by adding number of the tear variables is equal to the number of the re-
new nonlinear equations to the system. The nonlinear expres- maining nonlinear equations. It has been proved that these
sions for replacement are selected so that the net effect of each operations do not change the convergence rate of the Newton-
replacement is reduction of the number of nonlinear equations Raphson method. Numerical experiments have shown that they
in the system. The operations of this step are carried out in a may slightly improve the convergence rate of Broyden’s
Boolean matrix so that no formula manipulation is required. method.
2. Tearing of the linear subset of equations is carried out by The net effect of the tearing is reduction of the dimension of

Page 92 January, 1984 AlChE Journal (Vol. 30, No. 1)

the problem from n,where n is the total number of equations, tearing and as a result the computer time was reduced by
to the number of the nonlinear equations left in the system. 80%.
Since the computational effort is proportional to the square of We expect that the proposed tearing technique together with
the dimension in both the Newton-Raphson and Broyden’s the existing partitioning algorithms will play an important role
methods, quite substantial savings in computer time can be in making solution of large systems feasible and more effi-
achieved by using this tearing method. In one of the examples cient.
we have solved, the dimension was reduced from 39 to 17 by

INTRODUCTION ducible subsets. (“Irreducible subset” means a subset that cannot

be further partitioned.)
Steady-state simulation or design of chemical processes gives rise Detailed description of the principles of “partitioning” as well
to the need for solving large sets of nonlinear algebraic equations. as state of the art reviews of this field were recently published by
The need to solve such large systems became even more empha- Gustavson (1981),Sargent (1978), and Shacham et al. (1981).
sized lately, with the development of the equation-oriented flow- Irreducible subsets can be further decomposed by “tearing.”
sheeting or simulation programs (Shacham et al., 1981).Simulating Using this technique, the number of variables which should be it-
the operation of a plant using the equation-oriented approach, one erated on can be reduced. The subset of equations is rearranged
may need to solve a system containing tens of thousands of linear by tearing and then solved in one of the following forms:
and nonlinear equations. 1. The tear variables are guessed. New values for all (including
A system of nonlinear algebraic equation can be written: the tear) variables are calculated by solving each equation for its
output variable. (Output variable of an equation is a variable for
4x1 = 0 which the equation can be directly solved. When calculated, the
where f is the vector of functions and xis the vector of variables. output variable is a function of only the tear or previously calcu-
The Newton-Raphson (NR) method or one of its modificationsare lated variables.)
usually used to solve the system (Eq. 1). Two possible formulations 2. The tear variables are guessed. Some of the equations are
of the NR method are: solved directly for the new values of the nontear variables. The
residual of the rest of the equations is used to correct the value of
Xk+ 1 = X k - Jklf(Xk) (24 the tear variables.
and Tearing techniques have been used widely and successfully in
recycle calculations in the “sequential modular” process simulation
programs, but they have had only limited use in solution of systems
of nonlinear equations. There are two major difficultiesin applying
where k is the iteration number, fis the matrix of partial derivatives tearing techniques for nonlinear equations:
of f with respect to x (the Jacobian matrix) and pk = X k + 1 - Xk. First, tearing leads to some compounding of nonlinearities. The
Broyden’s quasi-Newton method (Broyden, 1965) is frequently compounded effect may deteriorate or even prevent numerical
applied for solving chemical engineering problems. This method convergence (Mah, 1972; Soylemez and Seider, 1973).
can be written as: The second difficulty is associated with output variable selection.
(3) Tearing requires that some or all of the equations be rewritten so
that one of the variables is an explicit function of the others. This
where Hk is the kth approximation to the inversed Jacobian matrix. cannot be always done. Even if it can be done it requires compli-
The Hmatrix is estimated for the first iteration, and after that it cated formula manipulations which are quite expensive in com-
is updated in every iteration according to Eq. 4. puting time (Hernandez and Sargent, 1979).
In this paper a new method for tearing systems of algebraic
(4) equations is introduced. This method does not require formula
manipulation, and it either improves or does not change the con-
where Yk = f ( x k + l ) - 4xk) vergence rate.
The direct application of Eqs. 2 or 3 for large systems is very In the next section the concept of tearing of linear subsets will
time-consuming and may be even impractical. The NR method be introduced. After that a method for eliminating nonlinearities
requires calculation of n2 partial derivatives (wheren is the number from some of the equations will be introduced. Several examples
of equations and unknowns) and storing and inverting the Jmatrix. will be given to demonstrate the advantages of the proposed
The use of Broyden’s method requires storage of the H matrix (n2 technique.
storage locations), and the update of this matrix requires about 3n2
multiplication operations (Dahlquist et al., 1965, p. 252). In addi-
tion, the initial Hmatrix is often estimated by calculating n2 partial
derivatives using numerical perturbation at the initial point and TEARING LINEAR SUBSETS
inverting the resultant Jacobian matrix.
Fortunately systems of equations that represent the operation Let us assume that the system of equations has been already
of a chemical plant tend to be sparse; thus, each equation will partitioned into irreducible subsets and the system (Eq. 1) repre-
usually involve only ten or less variables. This sparsity can be ex- sents one of these subsets. Let us further assume that this system
ploited for more efficient solution. contains both linear and nonlinear equations where the total
Large sparse systems can usually be “partitioned” into sets of number of equations is n, out of which m are linear equations. The
smaller problems which can then be solved sequentially, results system of equations can be rewritten:
from one problem providing data for the next. Since the storage
requirements and computation effort of both the NR and Broyden’s Ax= b (54
methods are proportional to n2, partitioning a system of equations
fl(X) =0 (5b)
into several problems with smaller n is clearly advantageous. It was
shown by Hernandez and Sargent (1979) that the computation cost where f i t R n - m is a subvector of f, bcRm is a vector of constants,
will usually be minimized by partitioning the system into irre- and A is an m X n matrix of coefficients.

AlChE Journal (Vol. 30, No. 1) January, 1984 Page 93

 ZU. The modified linear subset
1
1
-1

-2k5x5
Figure 1. Simple pipeline network and representative equations. Zb. Jacobian matrix of the unmodified system
r-2(kixl + kzxz + k 3 ~ 3+ k 4 ~ 4 )

Let us restructure the linear system (Eq. 5a) by reordering or by 2c. Jacobian matrix of the modified system
elimination so as to arrive at the following equivalent system: Flgure 2. Modlicatlonof the system of equationsof Figure 1 by tearing of linear
subset

where x1cRrnand x z ~ R " -are ~ subvectors of x, L is an m X m

lower triangular submatrix of A with nonzero main diagonal and When selecting the tear variables, one has to make sure that the
S is an mx(n - m )submatrix of A. resultant L matrix is not singular. Algorithm for selection of tear
If values are assigned to the components of the vector x2 (the variables is rectangular systems available in the literature (for ex-
vector of the "tear" variables), xl can be calculated from Eq. 6 by ample, Edie & Westerberg, 1971).
substitution since L is lower triangular. The solution of Eq. 6 for For complicated problems or for mixed systems of equations and
x1 can be formally written: procedures where analytical expressions for the partial derivatives
are not available,derivatives can be calculated by numerical per-
XI = L-l(b - SXZ) (7) turbation of the tear variables. Calculation of the derivatives this
The system of equations consisting of Eqs. 7 and 5b is equivalent
+
way requires n - m 1evaluationsof the system made up of Eqs.
7 and 5b. In the unmodified system (Eq. l),the vector of functions
to the original system (Eq. 5),but in the modified system only the
n - m tear variables are independent. Thus, the dimension of the
+
has to be evaluated n 1time, for calculating the derivatives nu-
merically.
problem has been essentially reduced from n to n - m. Slightly different formulations of Eqs. 7 and 10are obtained if
To use the NR method for solving this problem, the Jacobian a complete Gauss-Jordan elimination is carried out when trans-
matrix has to be evaluated. The Jacobian matrix for the original forming the matrix A into the matrices Land S.In this case L =
problem (Eq. 5) is the n X n matrix: Iwhere I i s the unity matrix and neither inversion of it nor mul-
r 1 tiplication by it is necessary. The computational effort associated
with the use of the two formulations is about the same, but this last
formulation is more convenient to use when the elements of JM
have to be expressed explicitely.
In the numerical examples that follow the Gauss-Jordan for-
mulation is used except in example 2. In all the examples, the
The (n - m)(n- m)Jacobian matrix of the modified problem, number of the tear variables was dictated by the number of non-
dflldxz,can be expressed using the chain rule: linear equations present. The tear variables were selected SO, that
the L matrix be nonsingular. Analytical derivatives are used in all
dfl - bfl I bfl dXl the examples and the inversed Jacobian at x = ~0 is used as initial
(9)
dxz bxz ~ X dx2 I estimate for Hoof the Broyden's method.
Substituting the value of x1 from Eq. 7 into Eq. 9 gives:
Example 1
Figure 1 shows a water distribution pipeline network and the
At this point the algorithm for using the NR method for a system system of equations that represents this network. (This example
containing linear and nonlinear equation that was modified by was taken from Ingels and Powels, 1964.) The parameters 91,
tearing can be summarized: 92. . .q5 are specified inlet and outlet flow rates, and the variables
1. Assume an initial value for the tear variables: (xz,~) x1 to 27 are the unknown internal flow rates. There are two laws
2. Compute x1.k from Eq. 7 and flxk) from Eq. 5b (where k is that govern flow rates in a pipeline network: 1.the algebraic sum
the iteration number; 0 at the beginning) of flows at each node must be zero; and 2. the algebraic sum of
3. Compute the Jacobian matrix using Eq. 10 pressure drops around any cyclic path must be zero. The applica-
4. Use NR iteration (Eq. 2a or 2b) to calculate the new a p tion of the first law leads to linear equations (the first five equations
proximation for the tear variables Xz,k+l in Figure 1);the second law gives nonlinear equations like the last
5. Check convergence. If converged, exit; otherwise, go back two equations in Figure 1. (For the sake of simplicity and clarity,
to step 2 a constant value for the friction factor was used in this ex-
It should be emphasized that the convergence rate of the NR ample.)
method is the same for the original and the modified system. The The variables x1 and zg were selected as tear variables, and using
reason for that is that the Jacobian matrix is not altered by tearing. Gauss-Jordan e l i a t i o n the linear subset of equations was brought
One can bring the original Jacobian matrix (Eq. 8)into the modi- into the form of Eq. 7. The modified form of this subset is shown
fied form (Eq. 10)simply by Gauss-Jordan elimination. (This was in Figure 2a. The Jacobian matrix for the unmodified system is
shown for example by Mah, 1972.)Thus, when started from the shown in Figure Zb, and that of the modified system is shown in
same initial point, the same intermediate values of Xk will be ob- Figure 2c.
tained in the original and modified problem. (Note that using the The two forms of this problem have been solved using the NR
same initial point means that ~0 has to satisfy Eq. 5a when the and Broyden's methods, starting with the initial estimate, Table
original formulation is used.) 1. The convergence of the solution methods for this problem is

Page 94 January, 1984 AlChE Journal (Vol. 30, No. 1)

 TABLE
1. INITIAL ESTIMATE
AND FINAL
RESULTSFOR PIPELINE TABLE
2. NUMBER OF ITERATIONS AND COMPUTATIONAL TIMES
NETWORKPROBLEM FOR EXAMPLE 2

Initial Estimate Solution Problem Solution No. of Computation*

Formulation Method Iteration Time(s)
21 690.0 515.24
22 696.3 871.06 Unmodified NR 3 1.578
23 280.4 455.16 Modified NR 3 0.313
24 551.4 245.63 Unmodified Broyden 5 3.043
x5 138.6 269.60 Modified Broyden 5 0.556
26 0.0 131.0
27 415.9 284.90 The computation time is based on a CDC Cyber 173 computer.

8 0 NR, modified

18

Figure 4. Pipeline network for example 2.

I I I I I I I
0 1 2 3 4 5 6 7
ITERATION NUMBER
Figure 3. Convergence of solution methods for example 1.

REDUCTION OF NONLINEAR EQUATIONS

The application of the tearing method, as presented previously,

shown in Figure 3.In this figure, the logarithm of the Euclidean is somewhat limited, since it requires some of the equations to be
norm of the residuals after each iteration is shown. completely linear. In this section we will show how the number of
The initial estimate used satisfies the linear subset of equations. linear equations can be increased and the number of nonlinear
As a result, the convergence of the NR method was identical for equations reduced in a system without changing the convergence
the modified and unmodified forms of the system. For Broyden’s rate of the NR method.
method there were slight differences in convergence for the two A nonlinear equation can be considered in general as summation
types of formulations. But this difference did not show any clear of three types of terms:
trend and it is too small to be of any significance. a. Constants
b. Linear terms: a constant multiplied by a variable
Example 2 c. Nonlinear terms: a constant multiplied by a nonlinear ex-
pression
A problem of larger dimension was selected to demonstrate the The following theorem will be useful in eliminating nonlinear
savings in computer time that can be achieved by tearing the linear expressions from equations:
subset.
This example is a pipeline network which is taken from Mah and Theorem 1
Shacham (1978).The network is shown in Figure 4. All the nodes
in this network are at the same elevation and all pipes are 30.4 m A nonlinear expression A x ) , which appears in one or more
long and 0.152 m in diameter. The equations representing this equations of the system can be replaced by a new variable: xn+l
network are similar to those of example 1, except that in this case and a new nonlinear equation: fn+ 1 = xn+ 1 - g(x) = 0 which is
the dependency of the friction factor on the flow rate was taken added to the system, without altering the convergence rate of the
into account. NR method for this system. (The proof for this theorem is given
Thirty nine equations are needed for representing this system, in Appendix A.)
out of which 22 are linear equations. This problem was solved in In chemical engineering problems there are many examples
both the modified and unmodified form using the NR and Broy- where the same nonlinear expression appears in several equations.
den’s methods. Initial estimates and final flow rates as well as ad- When carrying out steady-state material and energy balance on
ditional data for the problem can be found in the article by Mah a CSTR, the same nonlinear expression for the generation of say,
and Shacham (1978).Number of iterations and computation time species A, may appear in several material balance equations and
for the different methods are shown in Table 2. Convergence was in the energy balance equation. In such case this nonlinear ex-
achieved when the Euclidean norm of residuals became smaller pression can be replaced by a new variable. The total number of
than 10-4. equations will increase as a result of that; but if the number of linear
Comparison of the computational times shows that modification equations generated is larger than the number of nonlinear equa-
by tearing of the linear subsystem reduced the computationaleffort tions added to the system, the net effect is the reduction of the di-
to about l/5th of that of the original system for both the NR and mension of the problem. Also, even though the original system was
Broyden’s methods. Thus the great advantage of tearing the linear an irreducible subset, after modification it can often be fhrther
subset has been clearly demonstrated. partitioned.

AlChE Journal (Vol. 30, No. 1) January, 1984 Page 95

 f&) = x8 - l / x 7 =0

f9k) = X9 - %/x7 =0
Figure 7. Modifled form of the system of equations from Figure 5.
SELECTION OF NONLINEAR EXPRESSIONS FOR
REPLACEMENT

Replacing one or more nonlinear terms by a new variable can

be beneficial only if the number of equations that become linear not be any advantage in replacing this expression. The row and
by this replacement is higher than the number of new nonlinear column in which this expression appears can be eliminated from
equations generated. This simple rule can be used for selecting the matrix and the respective equation added to the subset of
nonlinear expressions for replacement. The selection can be done nonlinear equations.
by inspection only if the system is not too complicated. The selec- After the simplification, the search for expressions that can be
tion procedure can be, however, carried out by the computer. This replaced starts. At the first stage of the search, a single column is
analysis is carried out by the following method. considered at a time. The algorithm for this stage is the fol-
First the system of equations is scanned for nonlinear expressions, 1owing:
where each new nonlinear expression different from all the pre- 1. Select the next column with the highest frequency.
vious ones is assigned a symbol nel, nez . . . etc. The nonlinear ex- 2. Delete this column on a trial basis. If as a result the rank of
pressions can be compared, for example, by assigning a different more than one question becomes zero, the correspondingexpression
primary number to every variable and comparing the numerical should be replaced with a new variable. Delete this column, and
values of the expressions. Equal numerical value means that the any row whose rank becomes zero, from the matrix. Put the cor-
two expressions are equal. responding equations into the linear subset. Generate a new non-
After scanning the equations, an expression occurrence matrix linear equation and put it into the nonlinear subset. If less than two
is prepared so that it has the elements: equations have zero ranks as a result of the trial deletion, keep the
column in the matrix.
T = (til)=
I
1 if expression el appears in equation i
0, otherwise
The rest of the analysis is carried out by using only this Boolean
3. If there is more than one row of rank 1 and there are still rows
that have not been tested, go back to 1,or else start the second stage
of the search.
matrix, without the necessity to manipulate algebraic formulas. In many practical cases, this first stage of the search will result in
Expression frequency means the number of equations in which the elimination of all the columns and rows of the occurrence
an expression appears and equals the sum of elements in a column matrix. In other cases,higher order searches will have to be carried
of the matrix T. Equation ranks is the number of nonlinear ex- out.
pressions involved in an equation and equals the sum of the ele- The use of theorem 1with first-order search for simplifying and
ments in a row of T. solving a typical chemical engineering problem will be demon-
The Occurrencematrix can be simplified by using the following strated in example 3. Higher-order search method will be intro-
rules: duced after this example.
1. Any row of rank 0 can be eliminated and the respective
equation included in the subset of the linear equations, since this Example 3
equation does not contain any nonlinear expressions.
2. A column of frequency 1 indicates that the correspondent Figure 5 shows a system of nonlinear equations which was taken
nonlinear expression appears only in one equation, and there can by Carnahan et al. (1969). This system represents material and
energy balance in an adiabatic reactor for partial oxidation of
methane. In this system there is only one linear equation. The ex-
Equatl on pression occurrence matrix for this system of equations is shown
in Figure 6.

1
Equation 1 1 Using the two rules of simplification, the matrix of Figure 6 is
nunber 2 1
reduced to include only rows 1 , 2 , 3and 4 and expressions 1 and
3 1
2, Eq. 5 is placed in the linear subset and Eqs. 6 and 7 are placed
in the nonlinear subset.
4 1 1
Using the first-order search, two nonlinear expressions are re-
5
placed with new variablesx8 = 11x7, and xg = %6/%7,and two new
6 1 1
equations are generated: x8 - 11x7 = 0 and xg - xg/x7 = 0. The
1 1 1 last two columns of the Tmatrix are deleted so that the search for
Expression
frequency 2 3 1 1 1 1 replacement of the expressionsstops. Figure 7 shows the modified
form of the system of equations. It contains five linear and four
nonlinear equations. This system can, however, be further parti-
tioned, since x7 appears only in equations fs and fg, and 26 appears
only in fg. The first seven equations have to be solved simulta-
Flgure 6. Expression occurence matrix for the problem of Figure 5. neously, after Eq. 8 can be solved for x8 and Eq. 9 for xg.

Page 96 January, 1984 AlChE Journal (Vol. 30, No. 1)

TABLE3. INITIAL ESTIMATES
AND FINALRESULTS FOR EXAMPLE
3

Initial Estimate Solution~

0.208 0.32287
0.042 0.009224
0.048 0.046021
0.452 0.61817
0.250 0.0037165
0.340 0.57672
2.0 2.9779
0.50 0.33581
0.170 0.19367

‘I o NR 1 Equation

1 2 3 rank

Equation 1
nunter 2
1

1 1
1

2
1
3 1 1 1 3
4 1 1
5 I 1
6 1 1 2
Expression
frequency 3 4 2

ne -‘ X X 2/5‘ neZ’X2X3/S~ ne3’X3X6/S~

Figure 10. Expresslonoccurrencematrix for the problem of Figure 9.

ITERATION NUMBER new nonlinear equations according to the number of expressions

Flgure 8. Convergence of the solution methods for example 3. that have been replaced and put them in the nonlinear subset.
Return to the first-order search. If less than i f 1 equations have
zero ranks as a result of the trial deletion, keep this group of col-
The variables x8 and 29 were selected as tear variables, and the umns in the matrix.
system was decomposed by Gauss-Jordan elimination as in the 3. Delete the column with the lowest frequency from the search
previous examples. As a result of these operations, the problem’s group and replace it with the column of highest frequency that has
dimensionality was reduced from the original seven to two. not been included in the search group yet. Go back to 1. If there
The modified and unmodified forms of this problem have been are no more columns which have not been included already in the
solved using the NR and Broyden’s methods, starting from the +
search group, start the i first stage of the search.
initial estimate, Table 3. The convergence of the solution methods The search for replacement of expressions stops if one of the
for this problem is shown in Figure 8, and the values of the variables following conditions is satisfied:
at the solution are shown in Table 3. 1. All the columns have been eliminated from the matrix.
Figure 8 shows that the convergence is slightly better for the 2. The minimal rank is higher than the maximal frequency.
modified form when the NR method is used and better when the This algorithm for selecting nonlinear expressions for replace-
modified form of Broyden’s method is used. ment was tested with several practical problems. In all the cases,
Thus, tearing a system of nonlinear equations using the proposed the search stopped after the first or the second stage of the search
method may even improve the convergence instead of deterio- when all the columns have been eliminated from the occurrence
rating it. Such deterioration was observed when using other tearing matrix.
algorithms. The use of higher order search will be demonstrated with the
The change in the convergence rate of the NR method for the following example.
modified system is a result of partitioning, since this operation does
not preserve the convergence rate like the other operations do. Example 4
Figure 9 shows a system of nonlinear equations which represents
material balance on a chemical reactor (The Williams-Otto Pro-
HIGH-ORDER SEARCH FOR REPLACEMENT OF NONLINEAR cess). 91,41, . . .96 are specified inlet flow rates; k l , kz and k 3 are
EXPRESSIONS specified reaction rate coefficients; and X I , x 2 . . .xg are the un-
known outlet flow rates.
Sometimes the use of the simplification rules and first-order All the equations in this system are nonlinear. The expression
searches cannot eliminate all the columns from the expression oc- Occurrence matrix for the system is shown in Figure 10.
currence matrix. In such cases, the use of higher-order searches The matrix cannot be simplified since there are no rows of rank
(starting with second order) is necessary. The algorithm for an 0 or columns of frequency 1. Deleting a single column at a time
ith-order search is the following: zeroes only one equation rank; thus, first-order search cannot
1. Select i columns with the highest frequency. simplify the matrix. Deleting columns 1and 2 together zeroes the
2. Delete these columns on a trial basis. If, as a result, the rank ranks of Eqs. 1,2, and 4. These columns can be deleted and new
of more than i equations become zero, the corresponding expres- variables, x7 = x1x&2 and x8 = x 2 x 3 / S 2 , are defined together with
sions should be replaced with new variables. Delete this group of the corresponding nonlinear equations.
columns and any row whose rank became zero from the matrix. After this step, only column 3 and rows 3,5and 6 remain in the
Put the corresponding equations into the linear subset. Generate matrix. Deleting colum 3 zeroes the ranks of all the remaining rows.

AlChE Journal (Vol. 30, No. 1) January, 1984 Page 97

A new variable xg = xsx6/S2 can be defined and the search can be
stopped.
In the final form, the system of equations contains six linear and
three nonlinear equations.
-*
3x1
bgl
bx2
....-bgl
axn
1
7

NOTATION
dfidfa ....-df2 0
A = matrix of coefficients ax1 ax2 dxn
b = vector of constants
f = vector of functions
g = vector of nonlinear expressions
H = approximation to inversed Jacobian matrix
I = unity matrix bfnbfn ....-afn 0
J = matrix of partial derivatives (Jacobian matrix) ax1 ax2 dxn
k = iteration number dgzdgz ... -
agn
1
L = lower triangular matrix -ax1 ax2 axn - X=Xk

m = number of linear equations in an irreducible subset

n = total number of equations in an irreducible subset - - - -
ne = nonlinear expression
4 = inlet or outlet flowrate
P1,k gl(xk) + xn+l
T = expression occurrence matrix P2,k fdxk)
X = vector of variables
x1,xz = subvectors of x
= -

APPENDIX PROOF FOR THEOREM 1

Theorem 1. A nonlinear expression g(x) which appears in one Pn,k f n(xk

or more equations of the system can be replaced by a new variable,
x,+l, and a new nonlinear equation, f n + 1 = x,+ 1 - g ( x ) = 0, Pn+ -gZ(Xk)
which is added to the system, without altering the convergence rate - -C=X) - xn+ 1.

of the NR method for this system.

Proof. The elements of Pk as calculated from Eq. 2b are identical
for the original and the modified system. Assume, without loss of
generality, that one of the equations, say f l ( x ) , is a sum of two
+
terms, f 1(x) = g l ( x ) g2 (2). Equation 2b written for this system’
has the following form:

-
...-ag1 +-dg2
axn ax,
....-bf2
dxn

....-afn
ax, d x = q

(All

Page 98 January, 1984 AlChE Journal (Vol. 30, No. 1)

Mah, R. S. H., “Structural Decomposition in Chemical Engineering Springer-Verlag, Berlin (1978).
Computation,” 72nd AlChE Nat. Meeting, St. Louis (May 21-26, Shacham, M., S. Macchietto, L. F. Stutzman, and P. Babcock, “Equation
1972). Oriented Process Flowsheeting,” Comp. Chem. Eng., 6(2), 79 (1982).
Mah, R. S. H., and M. Shacham, “Pipeline Network Design and Synthesis,” Soylemez, S., and W. D. Seider, “ A New Technique for Precedence-Or-
Ado. Chem. Eng., 10,125 (1978). dering Chemical Process Equation Sets,” AIChE J., 19,934 (1973).
Sargent, R. W. H., “The Decomposition of Systemsof Procedures and Al-
gebraic Equations,” G. A. Watson, Ed., Numerical Analysis-proceed- Manuscript receioed April 2, 1982; reoisim received December 6, and accepted
ings, Biennial Conf., Dundee, Lecture Notes in Mathematics 630, 158, December 14,1982.

Fate of Solids Fed Pneumatically through a

Jet into a Fluidized Bed
Solid tracer particles were fed pneumatically through a jet into a fluidized bed IFlYENlA KECECIOGLU,
to simulate the feeding of solids via a pneumatic transport line into a fluidized-bed
reactor operating in the slugging-bed mode. The fluidized bed was defluidized WEN-CHING YANG, and
instantaneously at different times after the initiation of the tracer particle injection. D. L. KEAIRNS
The bed was then sampled layer by layer to provide the radial and axial concen-
tration profiles of the tracer. Regular and high-speed movies (1,000 frames per Westinghouse Research and Development
second) were taken to study the operation of the fluidized bed and the phenomena Center
of the gas-solid two-phase jet. Pittsburgh, PA 15235
Experimental results on solid mixing, jet constriction and slugging frequencies,
slugging bed height, slug length, jet penetration, and jet half-angle at three nominal
jet velocities of 52,37, and 25 m/s and corresponding solids loadings are presented.
Additional experimental results on jet constriction and slugging frequencies, and
slug volume (axial slug size) obtained for a wider range of jet velocities confirm
the hydrodynamic trends observed during the tracer particle injection experiments.
The results indicate that solids mixing increases, and well-mixed conditions are
reached earlier, with an increase in jet injection velocity. The obtained mixing
times were correlated successfully in terms of the excess gas velocity. The experi-
mental data on jet penetration and slug motion were satisfactorily correlated by
modified versions of existing theoretical relations. The modifications included
the effect of the injected solids on jet penetration and jet half-angle and also the
effect of our semicircular column geometry and single wallslug configuration on
the observed slug motion.

SCOPE
Fluidized-bed reactors have been widely acclaimed in the actor operating conditions, knowledge of the particle history
process industry for their advantageous characteristics of good inside the fluidized bed is essential.
mixing, high rates of heat and mass transfer, mechanical ro- The fate of solids fed pneumatically through a jet into a
bustness, and capability of continuous operation. Research fluidized bed depends on the particle size and density, the jet
seeking to harness these favorable reactor characteristics to the velocity, the solid loading in the jet, the bed height, the bed di-
efficient and environmentally acceptable processing of fossil ameter, the jet nozzle design, and the fluidization conditions
fuels (e.g., coal combustion, coal and oil shale gasification) has in the bed. Despite the importance of the fate of pneumatically
been undertaken intensively in recent years. fed solids to reactor performance, the complex phenomena af-
In gassolid fluidized beds, the gaseous reactant is introduced fecting solids mixing have not as yet been studied systemati-
uniformly through a distributor plate, while the solids are in- cally.
troduced via individual feed points. As a result, one has to rely A simple but effective technique was developed in this study
on the characteristics of the flow field to promote solids mixing to determine the extent of solids mixing and the deviation from
within the reactor. the ideal flow patterns in different regions within a fluidized-
A commonly employed method of feeding solids into a bed reactor by sampling the injected solid tracer particles upon
fluidized bed is pneumatic injection of the solids through a defluidization.
twephase gas-solid jet. If the gas-solid reaction is fast, unde- Experimental results on solids mixing and the related flow
sirable, localized, nonuniform temperature and concentration field phenomena in a slugging fluidized bed are presented at
conditions may appear around the solids feeding points in the three nominal jet velocities of 52,37, and 25 m/s. The experi-
fluidized bed. Such an occurrence may compromise the per- ments were performed at atmospheric pressure and ambient
formance of the fluidized-bed reactor. Therefore, from the temperature.
points of view of reactor design and selection of optimum re-

AlChE Journal (Vol. 30, No. 1) January, 1984 Page 99