Page 286 - Introduction to Computational Fluid Dynamics

P2: IWV
P1: IWV/KCX
0 521 85326 5
CB908/Date
0521853265c09
9.4 STONE’S METHOD
7. Go to step 1 with the next value of j (i.e., j = j + 1). May 11, 2005 15:41 265
8. Repeat steps 1–7 until j = JN − 2.
A similar procedure can be executed for sweep on the j direction. Finally, we
note that a procedure for simultaneous solution for three, four, or more consecutive
lines can also be devised but the associated algebra is very tedious. It will be realised
that if a simultaneous solution procedure is devised for all lines in a given direction,
one will have a procedure that is equivalent to the matrix inversion method for the
whole field.

9.4 Stone’s Method

As mentioned in Section 9.1, the convergence rate is sensitive to the structure of
the coefficient matrix. Stone [79] devised a whole-field procedure that reduces this
sensitivity. To apply the method, it is first necessary to change the 2D node address
(i, j) to the 1D address N. Thus,

N = i + ( j − 1) × IN, (9.30)

where N = 1,..., N max and N max = IN × JN. Equation 9.5 therefore can be
written as

AP N N = AE N N+1 + AW N N−1
+ AN N N+IN + AS N N−IN + Su N . (9.31)

In matrix form, this equation can be written as

| A || |=| Su |. (9.32)
Figure 9.1 shows the fully expanded form of Equation 9.32. Note that matrix A has
a maximum of five nonzero elements in each row.
The main idea in Stone’s method is to represent matrix A as a product of two
matrices, U and L. Thus, the L matrix (or lower matrix) is formed in such a way
that all entries above the diagonal are zero. The diagonal element is occupied by 1,
and positions of −AW and −AS in the A matrix are now taken by BW and BS
(say). Similarly, in the U matrix (or upper matrix) all elements below the diagonal
are set to zero; the diagonal elements are occupied by BP and elements occupying
positions −AE and −AN are replaced by BE and BN, respectively. The L and U
matrices are shown in Figure 9.2. Note that the size of L and U matrices is again
N max × N max .
Unfortunately, the product matrix |U|×|L| does not produce the A matrix
exactly. Instead, a matrix shown in Figure 9.3 is produced. This matrix has
two additional nonzero entries that occupy positions NW(i − 1, j + 1) and
SE(i + 1, j − 1). In terms of elements of the U and L matrices, the elements of