well distributed around the perimeter of the wiped wall, and the resulting irregular residence time distributions produce low Peclet numbers. The range given for one data point in Figure 4 is for five measurements a t the same conditions. Conclusions
Although the diffusion model may not be valid a t low values of the Peclet number, the experimental results show the magnitude of improvement which can be expected from a change in conditions. If the wiper speed is reduced, however, it may be necessary to increase the number of wipers to maintain complete lateral mixing of the film as well as a high rate of heat transfer. The effect of feed rate shows that, for a given concentration of the more volatile component in the effluent liquid, the feed rate can be increased with a less than proportional increase in the heat requirement because a lower Q will be required. The material balance equation is written for a constant diffusivity, but the question must be raised whether the diffusivity changes as the volumetric flow rate is diminished by vaporization. If the H and Pe functions are approximated by straight lines, H=mG+n Pe = p G + q then
D = - =UL -Pe
L2G
Pe H
L2
np
+ m q + (mpG + n q / G )
The value of D is therefore constant only when d - (nzpG dG
+ nq/G) = 0
that is, when
The curves in Figure 1 for intermediate values of the Peclet number are thus only approximate a t higher values of Q, and the use of a higher number of increments in 2 for the solutions of Equation 1 was not considered justified. This matter concerns only the intermediate Peclet numbers, since exact solutions for zero and infinite Peclet number are easily obtained. Nomenclature
a
= relative volatility
D
= diffusivity, sq. cm./sec. = length of wiped film evaporator, cm. = volumetric flow rate, cc./sec.
L G H
= holdup, CC. n, p , q =constants = Peclet number Pe = Peclet number for distillation, based on feed PeO velocity = fraction vaporized Q = second central moment of residence time dis892 tribution, secS2 = first moment of residence time distribution, sec. U = velocity, cm./sec. 5 = mole fraction of more volatile component = mole fraction in feed 20 2 = dimensionless axial length
m,
e
literature Cited
Fletcher, R., Powell, M. J. D., Computer J. 6, 163 (1963). van der Laan, E. Th., Chem. Eng. Sci. 7, 187 (1958). RECEIVED for review January 23, 1969 ACCEPTEDFebruary 28, 1969
I T E R A T I V E PROCEDURE FOR S O L U T I O N OF S Y S T E M S OF P A R A B O L I C AND E L L I P T I C EQUATIONS I N THREE D I M E N S I O N S H. G. WEINSTEIN, H. L. STONE, AND T. V. KWAN Esso Production Research Co., Houston, Tex. 77001
An iterative procedure has been developed to solve systems of three-dimensional elliptic differential equations. The two-equation algorithm has been tested on an oil reservoir flow procedure is fast, is insensitive to rounding errors, and in very difficult problems it yields a a procedure such as the alternating direction technique cannot. For the problem tested, of iteration parameters is easily calculated from the coefficient matrix.
MANY physical problems require the solution of parabolic
and elliptic partial differential equations. These equations are frequently nonlinear and involve complex geometries and boundary conditions. Because of these restrictions, analytic solutions of the differential equations are at present impossible. Instead, one must seek solutions of the finite difference approximations of the equations through iterative techniques.
and parabolic problem. The solution where a reliable set
Many iterative methods have been developed. Most of these, including relaxation (Southwell, 1946) and successive overrelaxation techniques (Young, 1950), require excessive computer effort because they converge rather slowly or fail to converge. The more implicit alternating direction iteration procedure (ADIP) (Peaceman and Rachford, 1955) converges faster than the relaxation and overrelaxation schemes and, in general, requires less computational work. VOL.
8
NO.
2
M A Y
1969
281
More recently a new iteration technique, the strongly implicit procedure, or SIP (Stone, 1968) has been developed. Stone demonstrated in the solution of a single two-dimensional equation that SIP achieved greater rates of convergence than ADIP on all problems tested except the simple model problem in which the coefficients in the difference equation were constant and isotropic. The advantage of S I P over ADIP appears to increase as the complexity of the problem increases. S I P was extended to the solution of systems of two-dimensional equations by Weinstein et al. (1969). In the present paper, the SIP algorithms for solution of single and multiple parabolic and elliptic equations in three dimensions are developed. A test on a two-equation problem is also discussed. Single Equation Case
For purposes of generality, a parabolic partial differential equation is chosen for investigation. The development for the corresponding elliptic problem is analogous, the only difference arising in the main diagonal terms of the resultant coefficient matrix. The equation to be solved is given in Equation 1:
Equation 1 is given by
1 - r i , j , k( P i , j , k - P i , j , k ) At
+
Qi,j,k
(2)
wherepi,i,k and P i , j , k are the values of the dependent variable a t the old and new time steps, respectively, and At is the time step size. Combining and ordering terms results in
+
zi,jtkPi,jtk-l
Di,j,kPi-l,j,k
Bi,j#i,j-l,k
Fi,j,kPt+l,j,k
+
Hi,j,kPi,j+l,k
+
+
Ei,j,kPi,j,k
Si,i,kPi,i,k+l= q i , j , k
(3)
where
(11 where x, y, and z are position coordinates, Q is a source term, and P is the dependent variable. K X , K Y , K Z , and are known coefficients which in general are functions of position and the dependent variable P . The solution of Equation 1 is sought through its finite difference representation. A grid network is defined as follows : i = 1 , 2 , ... I x i = iAx y j = jAy zk =
kAz
... J k = 1, 2, . . . K
(KZi,j,k+f
f
AxAy AxAy Az -- r i , j , k AZ At
etc., 1
j = 1, 2,
A typical grid element is shown in Figure 1. For the point (i,j , k ) the finite difference representation of
KZi,j,k-+)
1
r i p j v z i , j t k AxAYAz
The points (i,j , k ) of the grid are ordered in sequence such that i is swept through first (i = 1, 2, . . . , I),j second ( j = 1,2, . . .,J),and k last ( k = 1,2, , . , , K ) . If we order Equations 3 in the same sequence, then the over-all system of equations can be couched in matrix form-namely,
MP=i where
-
Pl,l,l P2,l,l
PI,l,l
,.
P1,Z.l
P= P2,2,1
P1,2,1
,PI,J,K_ Figure 1. 282
I&EC
Typical grid element
FUNDAMENTALS
6 is similar to P, with M
is shown in Figure 2.
Figure 2.
Figure 4. Matrix U
Matrix M
I and IJ between arrows indicate number of matrix elements between diagonals
Figure 5. Figure 3. Matrix I
If a direct solution of Equation 4 is sought by factorization of M into lower and upper triangular matrices L and U , it is found that I, has elements filling the entire region between the diagonals corresponding to 2 and E; and likewise, U is filled between the diagonals corresponding to E and S. For any practical values of I and J , this factorization would require excessive computer storage and calculation time. We therefore seek to modify M by N such that ( M N ) is factorable into L and U matrices which are sparse and of the form shown in Figures 3 and 4. The calculational and storage requirements of this factorization would be considerably less than those of the direct solution. With this modification, Equation 4 becomes
+
+
( M + N)P = tj+ Ni,
(5 1
Since ( M N ) is readily factored, Equation 5 can be solved easily if the right-hand side is known. Thus the following iteration scheme results: (N
+ N)P”+‘= 4+ Niim
(6 1
where m is the iterative level, To reduce rounding errors, a formulation suggested by Wilkinson (1963) and Peaceman (1969) is used. This is the
Matrix 1U
so-called “re_uidualform” and is obtained by adding and subtracting MPm on the right-hand side of Equation 6. The latter can then be put in the form:
(M
+ N)dP,+l
(7 1
=Rm
where the residual array Rmis defined by Rm=tj-,&m
and A
h
,jpm+l= p
h
+ l -
= change
pm
of P in the (m
+ 1)st step of iteration
In solving for the (m+ l ) s t iterate of @ (or d P ) , the right-hand side, involvjng only_Pm,:i known. When Equation 6 has con_verged,Pm+l= ,Pm = P; when Equation 7 has converged, am+’ = 0 and Pm = 6. The effects of the modification disappear, and the solution of Equation 4, which is the desired result, is obtained. Equation 6 (or 7 ) is the basic equation of SIP. The modified matrix (M N ) must satisfy the following relation : ( M + N ) = LU (8)
+
When the matrix multiplication is performed, a matrix of the form shown in Figure 5 results. N can be specified in many VOL.
8
NO, 2 M A Y 1 9 6 9
283
Combining Equations 7 and 8 leads to With vector $ defined as follows: e
Udpm+l=
v A
Equation 11 becomes
~ t=jR m For the point (i,j,k ) , Equations 13 and 12 are, respectively,
V i , j , k= d 2. , .3 h-1 (R 2., 3.,k m - a.1.3,k. V. . a~h-1-
- G,j,hVd-l,j,k)
(14)
- S i ,j,kdPi,j,k+lm+'
(15 )
bi,j,kVi,j-l,k dpi,j,km+'
= vi,j,k
- eiSj.kdpi+i,j,km+' f i,j , k d P i , j + l
,k"+'
-
Equations 10, 14, and 15 comprise the 3-D SIP algorithm. They form a determinate system for the solution vector dPrn+lwhen the following "boundary conditions'' are specified: zi , 3. , 1 - 8 t ., .3 , ~ - O i=l,2 I ; j = l , 2 ... J
...
O
i = 1, 2
, ,,
0
j = 1, 2
.. . J ;
Bi,l,k = Hi,J,k = Dl,j,k
= Fl,j,k =
I; k
=
1, 2
k = 1, 2
,,,
K
... K
The sequencing of points discussed above is not used on every step of the iteration. As for the two-dimensional case (Stone, 1968; Weinstein et al., 1969), it is alternated with another sequence in successive steps. During odd-numbered iterations the above sequence is used, and for even-numbered it,erations the sequencing is such that i is swept through first in increasing order (i = 1, 2, . . , , I ) , j second in decreasing order ( j = J , J - 1, . . ., 1), and finally k in decreasing order ( k = K , K - 1, . . , , 1). This last ordering creates nonzero elements in ( M N ) a t locations corresponding to the x-point,s and solid circle points instead of the open and solid circle points in Figure 1. The alternation of sequences tends to effect an over-all symmetry in (X+N ) in the two-step process and accelerates the convergence of the procedure. (Two additional independent orderings are possible, but those described above have thus far proved successful.)
+
Multiple Equation Case
The system to be studied is comprised of coupled, threedimensional parabolic equations. These equations are presented in Equation 16:
(16)
KX,p, KY,p, KZ,p, and r,p are known coefficients which, again, are generally functions of x, y, z , and the dependent variables P,, and Q, are source terms. Finally, n is the number of equations to be solved. After the terms in the finite difference representation of Equation 16 are combined and reordered, Equation 17 results:
284
l&EC
FUNDAMENTALS
It follows that the iteration procedure, the algorithm, and the forward and back solutions, Equations 9, 10, 14, and 15, respectively, all apply to the present problem when a, b, c, d, e, f, and g are replaced by n X n matrices and Ti, dP, and R by n-vectors.
where
Iteration Parameters
As before, the superbar indicates values a t the previous time step. Matrix notation simplifies Equation 17 as follows:
f
zi,j,kpi,j,k-l
Bi,j,kpi,j-l,k
+
Di,j,kpi-l,f,k
where for n
=
=[
Ei,j,kPi,j.k
H r.. i. , k p r., j .+ l , k f 2 , for example,
zi,jskaa
Zi,j,k
+
zi,j,kb”
1,
+
Fi,j,kPi+l,j.k
Si.j,kPi,j,k+l
1- amax =
+
= qi,j,k
=[
From a von Neumann error analysis (Stone, 1968) based on the single-equation model problem ( K X = K Y = K Z = constant, Ax = A y = Az = constant), one finds that small values of a attenuate the high frequency error components and large values of a attenuate the low frequency errors. For the most rapid convergence it is best to use a sequence of parameters in cyclical fashion. The convergence rate does not depend critically on the minimum parameter, so zero is arbitrarily chosen. For the it has been found best to use maximum parameter amax,
(18)
where
1,
p1=
Ei,j,kM. E i s j v k d
zi9j,kd
Ei,j,k
etc.
Ei,j,kh
zitj,kbb
Ei,j,kbb
-+-
KYAx2 KXAgZ
KZAx2 KXA2
KXAy2 KYAx2
ZAy2 +-KKYAz2
KXAz2 KZAx2
KYAz2 KZAy2
PZ = -
q i , j , k = (‘“ih)
p i , j , k= r i , j , k ) , pbi j ,k
p3=-+-
qbi,j,k
I
and a and b are symbols given to the two equations and unknowns. Equation 18 appears very similar to Equation 3. Indeed, the only difference is that scalar coefficients have been replaced by n X n matrices and the scalars P and q have been replaced by n-vectors. Because of this similarity, the extension of the single equation S I P to t,he multiple equation case is straightforward. The grid points and corresponding difference equations (18) are sequenced as first discussed above. The over-all system of mabrix equations can again be placed in the matri; form of Equat’ion 4. Now, however, the (i,,i,k ) element of P is the n-vector P i , { , k ;for 6 it is q i , j , k , and the elements of 111 become the n X n matrices which were illustrated for n = 2 following Equation 18. With the present problem assuming the form of Equation 4, it follows that the basic equation of the 3-D multiple equation S I P will be the same as that of the single equation procedurenamely, Equation 6 (or Equation 7). The modifying matrix N must again be defined such that ( M N ) is factorable into sparse matrices which have the configurations shown in Figures 3 and 4. Also, matrix multiplication of L and U yields a matrix of the form shown in Figure 5. Whereas in the single equation case the elements of L, U, and ( M N ) were all scalars, t’hey are now n X n matrices for the n-equation case. Defining N analogous t o the scheme described in the previous section results in transforming C, G, A , W , T, U,and a p ( p = 1 , 2 , 3 ) into n X n matrices, where for n = 2, for example, a p is given by
+
+
0
ffbp
p
is a transmissibility ratio. The maximum of p over the grid
is given the symbol pmax. For 3-D problems,
max
Pmax =
[ p i , PZ, ~ 3 1
(over grid)
I n evaluating pmax and the minimum in Equation 19, extremes of pp = 0 or pp = m ( p = 1, 2, 3 ) are avoided. Equation 19 is the same as that used to calculate the minimum ADIP parameter for a homogeneous single-equation problem with rectangular grid elements. [This minimum parameter is derived by performing a von Neumann error analysis on the normalized ADIP equations. The normalization is indicated in Appendix B of Douglas et al. (1959) and the error analysis is discussed by Varga (1962).] For cases where p varies throughout the grid, local values of omaxare calculated and the maximum over the grid chosen. Generally, 4 to 10 parameters are used in a cycle. They are arranged geometrically between minimum and maximum as follows: (I - a,) = (1 - amax)m’(C-l)m = 0 , 1, . . , , (C - 1)
(20 ) where C is the number of parameters in a cycle. The parameter sequence given by Equations 19 and 20, although derived for the single-equation model ADIP problem, has proved effective in the solution by S I P of multiple-equation three-dimensional reservoir problems, the same parameters being used for all equations and planes. If the sequence causes divergence of the iteration, (1 - amax) as given by Equation 19 should be multiplied by a factor of 2 to 10. On the other hand, if the iteration is converging, but slowly, (1 - a m a x ) should be divided by a factor of 2 to 10. VOL.
8
NO. 2 M A Y 1 9 6 9
285
Evaluation of SIP
decrease if the tolerances were loosened, but with looser tolerances the solution would become less reliable.) Thus, the fully double precision ADIP program required over twice as much computer time as the single precision S I P program. The failure of single precision ADIP to converge on a problem with pmaxas large as 645, which was the value here, appears to be largely the result of rounding error. This is evidenced by the initial sharp drop of residuals in a time step, followed by essentially no change in their magnitudes. Double precision has alleviated some of these difficulties, as indicated above. It appears that SIP is not as critically affected by rounding error as ADIP. It is seen in Table I1 that the iteration time (milliseconds per iteration per block) is larger for SIP than for ADIP. However, because SIP, in general, can converge in fewer iterations, the latter requires less computational effort than ADIP for solving the same problem.
Extensive testing has been carried out on the two-dimensional SIP procedure for both single and multiple equation problems. These tests are reported by Stone (1968) and Weinstein et al. (1969). Based on the encouraging results of the tests, the three-dimensional S I P development was undertaken. The original version of S I P in three dimensions was derived to solve systems of equations. The one-equation algorithm has not been tested, but as a result of the two-equation test reported here, it is believed that SIP in three dimensions is just as effective as SIP in two dimensions. Using the IBhI 360/65 computer, the test was performed on a reservoir flow problem in which water was displacing oil from an oil reservoir. The two equations to be solved were material balances on oil and water, both of the form given in Equation 16. The data on the reservoir and fluids are presented in Table I. The problem was solved using both the SIP and ADIP iterative procedures, and results of the two methods are compared in Table 11. In a single precision program, ADIP failed to converge to a solution. When a double precision ADIP algorithm was embedded in a single precision program, convergence was achieved, but the desired tolerance could not be met on the first time step. The L1 norm (sum of absolute values) of the oil residual could be reduced only SO%, not the specified 90%. Only by running in a fully double precision program could ADIP meet the convergence criterion on all time steps. On the other hand, SIP in single precision met the tolerance on all time steps with no difficulty. Optimum parameter sequences were determined for ADIP, while the S I P parameters were chosen by Equations 19 and 20. A work ratio, defined as the ratio of ADIP computer time to SIP computer time, was calculated in order to compare the two iterative procedures. The work ratio calculated for the present problem was 2.20. (The work ratio would tend to Table I.
Conclusions
The strongly implicit iteration procedure has been developed to solve single equations and systems of equations in three dimensions. SIP has been found to have the following characteristics: 1. On the problem tested, SIP required significantly less computational effort than ADIP. 2. SIP is less subject to rounding errors than ADIP. 3. The iteration time of SIP is modestly larger than for -4DIP; however, SIP generally requires fewer iterations to converge. 4. For the problem tested, Equations 19 and 20 provide a reliable sequence of iteration parameters.
In conclusion, one can state that SIP is an efficient numerical iterative procedure for solving parabolic or elliptic equations in three dimensions. Appendix.
Physical Characteristics of Water-Oil Reservoir Problem
Reservoir length, ft. Reservoir width, ft. Reservoir height, ft. Grid blocks in s-direction Grid blocks in y-direction Grid blocks in z-direction Porosity, % Oil density, g./cc. Water density, g./cc. Oil viscosity, cp. Water viscosity, cp. s-Transmissibility, darcy-ft. y-Transmissibility, darcy-ft. ,+Transmissibility, darcy-ft. Pmax
Oil production rate, STB/day Water injection rate, STB/day Simulation, months
Derivation of 3-D Algorithm
A study of Equations 8 and 9 and Figures 3,4, and 5 , leads us to Equation A-1 :
1320 1320 13.5 10 10 7 6.7-21 0.81 1.0 1.1 0.62 1.7 X lod6 - 4.3 X lo-* 1.7 X 10-6 - 4.3 x 10-2 1.5 X 10-3 - 28 645 14.7 (i = 10,j = 10, k=1-7) 15 (i = 1, j = 1, k = 1 - 7) 0-62
1.
2.
3. 4. 5. 6.
7.
Table II. Three-Dimensional, Two-Equation Problem
Iteration Procedure SIP (S.P.)“
ADIP (F.D.P.)”
Parameter Max SIP (Min A D I P )
0.9999 2 x 10-2 (1st T.S. 10-4)
ReducSimulation in tion Period, Months
L.1 Norm
90% 90%
No. of
No. of
Time Steps
Iterations
Av. Itns./Time Step
Itn. Time, Maec./Itn./ Block
Total Time, Min.
Work Ratio
33 31
167 470
5.1 15.2
4.75 3.71
9.25 20.36
2.20
62 62
S.P.indicates “single precision” ; F.D.P., “fully double precision.” 286
ILEC FUNDAMENTALS
N P
8.
modifying matrix dependent variable vector of dependent variables right-hand side of difference equation, defined after Equation (3) vector of right-hand sides source term residual vector of residuals stock tank barrels time upper triangular factor of ( M N ) Cartesian position coordinates iteration parameters increments in time increments in position coordinates coefficient in differential equation transmissibility ratio
9
9.
4
10.
6
11.
Q
12.
R R STB
13.
The elements of & a i , j , k , h i , j * k , c i , j , k , di,j,k-and U-ei,,,k, gi,j,k-lie on the diagonals displayed in Figures 3 and 4. As an example of the method of solution, the element a i , j , k is chosen. Substituting Equations A-1-2 and A-1-12 into Equation A-1-1 results in
fi,j,k,
a.1 . 3‘, k =
Zi,j,k
- ~ ~ a i , j , k e i , j . k - -l Qsai,j,kfi,j,k-I
t
+
U
Y,
5,
Ql, 0 2 ,
ff3
At AX> AY,
r
P
or ai,j,k
=
1
+
SUBSCRIPTS
Zi,j,k ~zei,j,k-i
+
013 f i , j + 1
equation and variable designators position increments iteration level
P, a, h i, j , k
ff,
The remaining elements of L and U are obtained in similar fashion. Equation 10 presents the final algorithm.
m
SUPERSCRIPTS Acknowledgment
The authors are indebted to many individuals a t Esso Production Research Co. for their efforts in the development of the three-dimensional SIP. Particular thanks are given to B. L. Buzbee and H. J. Welge for their help in implementing and testing the procedure.
a,P, a, h m
-
equation and variable designators
= iteration level - old time level
literature Cited
Douglas, J., Peaceman, D. W., Rachford, H. H., Jr., Petrol. Trans. 216, 297 (1959).
Peaceman, D. W., private commiinication, 1969. Peaceman, D. W., Rachford, H. H., Jr., J . SOC.Znd. A p p l . Math.
Nomenclature
a, h, c , d , e, f, 9
A , c, G, T, U, W B, D,E , F , G, S , Z
C
e . dP
= elements of factors of
(M
+N )
= coefficients of modifying matrix = coefficients in difference equation = number of parameters in a cycle = change of P over an iteration = vector of dP
= number of grid elements in 2, y, z directions K X , K Y , K Z = coefficients in differential equation L = lower triangular factor of ( M N) M = coefficient matrix n = number of equations
1, J , K
+
3, 28 (1955).
Southwell, R. V., “Relaxation Methods in Theoretical Physics,” Clarendon Press. Oxford. 1946. Stone, H. L., SIA’M J . iumerical Anal. 6, No. 3 (September 1R6A)
Vaiga:’R. S., “Matrix Iterative Analysis,” Prentice-Hall, Englewood Cliffs, N. J., 1962. Weinstein, H. G., Stone, H. L., Kwan, T. V., “An Efficient Iterative Procedure for Solution of Two- and Three-phase Reservoir Flow Problems,” unpublished manuscript, 1969. Wilkinson, J. H., “Rounding Errors in Algebraic Processes,” Prentice-Hall, Englewood Cliffs, N. J., 1963. Young, D. M., “Iterative Methods for Solving Partial Differential Equations of Elliptic Type,” doctoral thesis, Harvard University, 1950. RECEIVED for review February 5, 1969 ACCEPTEDFebruary 5, 1969
VOL.
0
NO.
2
MAY
1969
207
