- Email: [email protected]
Comparison of two splitting algorithms for 2-D modelling of advection-dominated solute transport in groundwater
Comparison of two splitting algorithms for 2-D modelling of advection-dominated solute transport in groundwater Gnanadesikan Ramanujam and Antonis D. Koussis Department of Civil & Enviromnental Enghwering, Vanderbilt Universit.v, Nashville, Tennessee 37235, USA This paper reports on two alternating directions splitting algorithms, one cast in the Local One Dimensional (LOD) and the other in the standard Fractional Time Stepping (FTS) format. Both algorithms employ an adaptation of the pseudo-viscosity technique and use the principal directions of transport as coordinates. In the LOD approach, only one dimension appears in each splitting step. In the FTS method, computation is implicit in only one dimension but other dimensions are present in each splitting step. Stability and accuracy constraints on the grid parameters are diagrammed for both schemes. Test problems in one and two dimensions are solved. The analysis shows that the slight accuracy advantage of the FTS format over the LOD format must be traded off against tighter grid design prescriptions and, frequently, also against higher computing cost. Key Words:
Numerical, computation, transport, solute, advection, dispersion, model, groundwater, fractional step, splitting, alternating directions.
INTRODUCTION Fractional steps or splitting methods I were devised to increase efficiency in the numerical solution of multidimensional problems. By treating each spatial dimension implicitly in separate sub-steps, the single, large system of equations of fully implicit methods is decoupled and a series of smaller systems with matrices of smaller bandwidths results that are simpler to solve. In the Peaceman-Rachford Alternating-Direction Implicit (ADI) scheme ~ this is achieved by representing terms in a single dimension implicitly; terms in the remaining dimensions are approximated explicitly. Typically, the solution is advanced over one time step in sub-steps according to the dimensionality of the problem (6t/2 in 2-D, &/3 in 3-D). Yanenko ~ has also proposed a specialized ADI algorithm, which he terms a splitting scheme, that decomposes the governing equation to a number of strictly one-dimensional equations. Lapidus and Pinder-' refer to it as the Local One Dimensional (LOD) approach. We report on the comparison between two related ADI algorithms, one cast in the LOD form and the other in the traditional ADI form in which the space difference operator is multi-dimensional at each fractional step. The LOD scheme is derived in Yanenko 1 and Lapidus and Pinder ~ as a fractional steps method, but it can also be viewed as an approximation in two spatial process-steps, each operating over the entire time interval. Henceforth, we will refer to the standard ADI algorithm as Fractional Time Stepping or FTS method in distinction to the LOD methodology. Both algorithms employ an adaptation of Paper accepted January 1991. Discussion closes February 1992.
© 1991 Elsevier Science Publishers Ltd.
the 'pseudo-viscosity' technique 3, due to Cunge 4, which we term the Matched Artificial Dispersivity (MAD) scheme. Having developed a 2-D LOD groundwater solute transport model 5, we now explore whether retaining both spatial dimensions in each sub-step improves accuracy, and if so, at what cost and under what conditions. The algorithms are designed to integrate the mass balance equation along the principal directions of transport 6. These are the streamlines and a set of corresponding orthogonal coordinates (in media with isotropic hydraulic conductivity, the equipotentials). In this formulation, the transport equation has a single advective term and no cross-derivative terms. There is thus, principally, no difference in computing on Cartesian or on curvi-linear grids. Additionally, transverse numerical diffusion due to non-alignment of grid and flow lines is eliminated. Because the cost of generating the grid nodes on the curvi-linear coordinate system is substantial, the principal directions approach is limited in practice to steady flow fields. It is further assumed that the coordinates' curvature is small. Our interest here is in solute transport in groundwater, where the assumption of steady state flow is common practice 7. G O V E R N I N G EQUATION In many applications, transport of solutes in groundwater can be modelled adequately in two dimensions. With the principal directions of transport as coordinates, the mass balance equation for a conservative dissolved chemical
Adv. Water Resources, 1991, Vol. 14, No. 4
183
Comparison o f two splitting algorithms f o r 2-D modeling. G. Ramanujam and A. D. Koussis
species reads s ?C ?t + u?C.?x = u[?(x,?C/?x)/?x + ((:~.?C.?y)/(y] (1
in which C is the solute concentration, u is the linear pore velocity, x is the curvi-linear coordinate in the direction of flow, y is the coordinate transverse to the flow direction, :z., and :ty are the respective longitudinal and transverse dispersivities, and t is time. If the dispersivities vary so slowly in space (s) that terms (?:t//Ps) (?C/?s) are small compared to Di~ 2 C/~s 2, where D/= ucq, (1) takes the form £C/(t + u?C/?x = Dx 72C/?x2 -}- DyO2C/?y 2
(2)
The longitudinal and transverse dispersion coefficients, D., and D~., are not constant as u = u(x, y), in general. Dispersivities, however, are usually calibrated as constants (as in Ref. 9) because of uncertainty about the field conditions (geometry, flow field, source and observations). Furthermore, to describe the transport of solutes undergoing linear, reversible, equilibrium absorption, u is divided by a 'retardation factor 's. The numerical developments are based on equation (2) with constant coefficients. Space variable u, Dx and D~. can be accommodated by local averaging over grid elements. However, constant grid spacings are implied in all references to truncation errors.
advection dominates, the solution of the pure advection equation is a first approximation to the solution of the complete parabolic equation in the flow direction. On this premise, Cunge a develops an explicit, two-point finite difference (FD) approximation of the advection equation by modelling advection by a time-centred, first order-inspace accurate FD approximation, and by off-centring the time derivative in space (weighting coefficient 0 < 0.5). The discrete model of (3a), with L,. = - ua.,, thus reads ( o / , s t ) ( c * _ ,.~ - c 7 - ~ . j )
+ [(1 - o ) / ~ t ] ( c * j
+ (, 2`5x) [(c,*4 - c ~ _ ~ . j ) + ( c * j - c *
- c,*j) ~4)] = 0
where 3x is the longitudinal increment, ,5t is the time increment, and i, j are space indices, and * indicates the preliminary solution after one ,5t due to longitudinal transport. Rearrangement yields the following directionally explicit FD scheme, in which C = u,st/,sx is the Courant number, C*.=
K , C * _ ~ . j + K , _C ~ _ t . j
+ K3C*
Following the basic concept of the LOD method, equation (2) is split into two one-dimensional equations, here over one time step: (C* - C")/bt = LxC';
(3a)
and
K~ = (C + 20)/(2 + C - 20)
(5a)
Kz = (C - 20)/(2 + C - 20)
(5b) (5c)
The truncation error of this approximation gives rise to numerical diffusion. Cunge 'tunes' this numerical diffusion, to account for the initially neglected longitudinal dispersion term, and obtains second order accuracy in space. He accomplishes this by matching the numerical diffusion coefficient to the dispersion coefficient, the former being derived from a Taylor series expansion to second order. The match specifies the value of the spatial weighting coefficient 0 in terms of the grid Peclet number, P = u6x/D~:
( C n + 1 _ C * ) / 6 t ~- L y C ~
(3b)
where L x = - u6 x + Dx6 z, L~. = Oy62 and ,5i is a first order and ,52 is a second order difference operator. C* is a preliminary solution at the end of one time step ,5t due to longitudinal transport alone. This result is then used as initial condition for the computation of lateral dispersion, again over one 6t. The spatial derivatives are evaluated at some time level n ~< ~ ~
184
(5)
where
K = (2 - C - 20)/(2 + C - 20)
T H E LOD S C H E M E
(4)
Adc. Water Resources, 1991, Vol. 14, No. 4
0 = 0.5 - P - 1
(6)
The scheme is unconditionally stable for physically realizable processes, i.e. those with positive diffusion/ dispersion ~t. The Courant and grid Peclet numbers are subject to mild constraints that derive from the requirement that the scheme's coefficients be positive. Unphysical oscillations are thus eliminated and physically reasonable results ensured'*'5; this control of numerical dispersion does not interfere with the scheme's ability to compute steep fronts. Because of the time-centring of advection, and hence of dispersion (interpreting C* as C "+~ from longitudinal transport, or as C "+~/" in the FTS format), the scheme is also second order accurate in time. In the transverse step, second order accuracy is obtained by centring the second derivative in space and time 1, i.e. (C~.~ ~ - C*j)/6t = Dy(`5~I/.;C "+1 + 6~1,,~C*)/2
(7)
in which 62C ] 14 = (Ci.j+ t - 2C/,j + C~.j_ 1 )//~j2
(8)
Comparison o f two .splitting algorithms j b r 2-D modeling: G. Ramanujam and A. D. Koussis
and 6y is the transverse space increment. As noted, the analysis of truncation error is simplified when starred values are interpreted as concentrations at t = {n + 1/2),3t. THE FRACTIONAL TIME STEPPING ALGORITHM To select an FTS algorithm in which to incorporate the MAD scheme, we considered a number of alternatives presented by Frind 1-' and by Burnett and Frind t3. We decided on the Peaceman-Rachford scheme t't'~ because of its second order accuracy and simple structure. This scheme consists of two steps Step 1
(C,+ 12 _ C")'(5t/2) = LxC "+ 1;_,+ L~.C"
(9)
Step 2(C "+1 - C " + I ~ 2 ) ( 6 t / 2 ) = L x C " + I ! Z + L ~ . C "+~
(10)
We modify this scheme to employ a special form of the MAD scheme in the first half-step in which advection is evaluated at t "+ ~ia. And in further contrast to the LOD method, we retain and evaluate transverse dispersion at t". Evaluation of the space derivatives is thus centred in time over the first half-step (at t "+ t4), which gives second order time accuracy. The temporal derivative is replaced by a weighted sum (weighting coefficient 0) of two FD approximations of the time derivative over &/2 at nodes i - l , j and i,j. The second sub-step is also centred in time, at /,,+3 4. The scheme's equations are Step 1: L_
l/z = lr. ~ , , ~/-r~,i - +1,j
, 1,j) Ci-
+(1 - O ~, (, _C, ,"+ j i'e
-
(il)
cT, j)]/(&12)
Step 2: L ~ C "+ "'- + L y e "+~ = C,")~ - C','.~:
1/2,/,(}
'9
~, t t l - )
(12)
where L.~=-uax:
for step l,
Ly=Dy6~
(13a)
and Lx= -uf.~+Dx6j;
Ly=Dy6~
for step 2. (13b)
Clearly, the solution at each fractional step involves more effort than in the LOD scheme because the difference equations contain terms from derivative approximations in the x- and )'-directions. But for the same reasons the FTS scheme is expected to be more accurate. After collecting terms, we derive for step 1 C~' ~ x,2
I
n
..t
n
t
= K , C , _ t.j + K;C'~+~.') + K3Ci.j -f- KaC'/.j+ ~ + K'sC'~.~-,
(14)
where K't = 2 0 / [ 2 ( 1 - 0) + C ]
(14a)
K2 = (C - 2 0 ) / [ 2 ( 1 - 0) + C ]
(14b)
K ; = [2(1 - 0 ) -
(14c)
and (14e)
& = D~.&;(dy):
Because of the implementation of MAD, the solution of the first step in this formulation is explicit, and thus simple. In deriving equation (14), longitudinal dispersion was neglected. Yet, upon expanding the grid functions in equation (14) in Taylor series around (xi, yj, t "+ ~s-,) to second order, we obtain a transport equation with the numerical longitudinal diffusion coefficient D v = udx/2 - udxO + ua~SU4
To account for longitudinal dispersion, the numerical diffusion coefficient is, again, matched with the physical dispersion coefficient, D., = D s. As before, the match is achieved through the selection of 0; however, in contrast to the LOD match condition (6), in the FTS method 0 depends on P and on C, 0=0.5-
P-~ + 0.25C
K;, = K ; = i , , , / [ 2 ( 1 - 0) + C ]
(14d)
(16)
As in the L O D scheme, only two points are employed implicitly in the longitudinal direction in the first step of the FTS method. Since the second derivative is modelled as an adjusted truncation error of the advection equation, no downstream boundary condition is imposed, nor can one be imposed. This is of little concern in advectiondominated cases, but limits the method's applicability when diffusion is significant. Moreover, it has the advantage that computing time can be economized by limiting the domain longitudinally to the effective transport area at any instant. In the second step, equation (10), the dispersive terms are approximated by the second-order-accurate, spacecentred F D ratios of equation (8). In attempting to resolve the 'smearing' versus overshoot dilemma arising in advection-dominated transport conditions ~4, we experimented with alternative advection schemes. Some were discarded due to deficiencies in accuracy, others because they imposed too restrictive stability constraints (e.g., Leonard's is multi-point, asymmetric Q U I C K scheme). We decided on a 'diffusion-corrected', two-point upwind advection scheme, which we implemented as follows. Upon substituting in (10) the two-point upwind scheme for the advection and space-centred differences for the dispersion terms, the grid functions are expanded about (xi, yj, t "+ k,2) in Taylor series to second order. We obtain again an advection-dispersion equation, except that the effective longitudinal dispersion coefficient is D.~+O.5uSx-O.25u26t. This distortion is rectified by adding to Dx the terms - (ucSx2) + (u2ijt/4). The corrected coefficient of longitudinal dispersion, (Dx).... = D ~ (uax/2) + (u2Dt/4), nullifies the numerical diffusion of the first order accurate approximation for advection, The second-order accurate scheme for the second step is then: ,+1 1 + (2 + 2/ly)C~,;: 1 - ~'~.'-i,j+ ,, r , + l z -l~yCi.~t-,n+
& . l / J 2 ( 1 - 0) + C ]
(15)
1,2
= (2 - C - 2 # . < ) C ' ~ ~"-' + s i . , ~ + 1,j
+ (C . . .l.l x J ( ~
+ lie
i - 1, j
Adv. Water Resources, 1991, Vol. 14, No. 4
(17)
185
Comparison o f two splitting algorithms Jbr 2-D modeling. G. Ramanujam and ,4. D. Koussis
where :tx = (DO ..... 6t/dx 2 = (D x - 0.5udx + 0.25ue 6t)&/6x 2
= ( C P ) - 0.5C + 0.25C 2
(17a)
and as before, ll,. = Oy(t/~) '2, (14e). Because in (1"}) only the transverse dispersion term is implicit, the problem is effectively one-dimensional. The coefficient matrix of the system of linear equations is tri-diagonal and can be solved efficiently by the Thomas algorithm. Accuracy and Stability
In the previous section, both schemes have been shown to be second order accurate in space and time. This has been ensured by time-centring, by lonigitudinal dispersion matching or correction, and by direct second order approximation of the transverse dispersion terms. This section addresses grid restrictions dictated by stability and by physical feasibility of solution. Physically feasible behaviour imposes positivity on the schemes' coefficients, to ensure that positive inputs cause positive responses. This is necessary for the computationally explicit sub-steps; in fact, conditions more restrictive than those imposed by stability may result. The constraints are summarized in diagrams delineating grid parameter regions that ensure meaningful results, which we term feasibility diagrams and feasibility regions, respectively. Both steps of the LOD scheme are unconditionally stable. The stability of the first step is ensured simply by D~. > 0, or 0 < 0.5, which is automatically satisfied for Dx > 0. And the second step has the unconditionally stable Crank-Nicolson format. Only a physical feasibility analysis is therefore required. Restrictions on the feasible region are imposed solely by the first step and from the requirements that K., > 0 and K 3 > 0 in equation (6) (K ~ > 0 always) 4. They are displayed in Fig. 1, along with the stability restriction 0 ~ 0.5 and the limitation imposed by advection-dominated behaviour P/> 2. In the FTS scheme both steps impose limits, Physical feasibility dictates constraints for step 1, numerical
C
C :
2 :~/, ~p_2~i.
0~<0.5+C4
2-20
=o
'
t N~
DN
> o - 05
0
(18)
C~<40
(19)
For equation (14) to give meaningful results, K'~, K ~ . . . K'5 should not be negative. Substitution of equation (16) for 0 in the denominator of these coefficients yields i-1+(2/P)+(C/2)], which is always positive. The numerators are analyzed to determine the sign of K'~, K; • " K'5. K'~ is positive because the minimum value of 0 is C/4 (when P = 2) and C > 0. K~. and K'5 are positive because D~. is positive. To show that K', should not be negative, we consider a zero initial condition and a positive input to the system. Since the downstream concentration response must be positive and only the second term on the right hand side of equation (14) remains, K"2 must be positive, whence c > 20.
(20)
or
0~1
(21)
The conditions (19)--(21) are most restrictive, and thus define the feasibility region for the first half-step shown in Fig. 2. The second step of the FTS method is conditionally stable. When the stability of equation (17) is analyzed by the von Neumann method, the amplification factor 3 obtained is G = [2 - C - 2px + (C + 2p x) cos(p6x)
05
Fi 9. 1. Grid parameter feasibifity diagram f o r the L O D scheme
186
C1>40-2
In fact with the equal sign, equation (18) applies to pure advection, P--+ .m, the correct solution of which is obtained for C = 2, 0 = 1. The scheme is thus stable for any P. That transverse dispersion enhances stability follows from the amplification factor of the yon Neumann stability analysis, which is maximized when Dy is eliminated (by choosing the y-spatial frequency that minimizes the denominator). The feasibility constraints for step 1 are obtained by examining the coefficients in equation (14) and the matching condition (16). Experience with the MAD method ~6"~7 suggests that P = 2 is a meaningful lower limit to characterize advection domination. Application of P > 2 to (16) yields
2(-0)>0, femible
or
For a negative step input from 1 to 0, after the input has ceased only the third, fourth and fifth terms on the right side of equation (14) remain, the last two of which are always positive. On a grid with a sufficiently high transverse resolution we may assume that C~'+~.j~ C7+ ~,j+ ~ ~ C'/+ ~.j_ ~ = 1, which effectively eliminates + 1/9 transverse dispersion. We then have for C ~ l~+ ~.'f > 0
ib [ ,
stability for step 2. In step 1, stability is satisfied automatically through the matching D x = D x, but 0=0(C, P). This can be seen from the positivity of physical dispersion coefficients, D., = D.v > 0, which can be interpreted as a stability constraint~: with equation (15), it yields
Adv. Water Resources, 1991, Vol. 14, No. 4
- \ f ~ - 1C sin(p6x)]/[2 + 2lty(l - cos(q6y)] (22)
Comparison
of two splitting algorithms jor 2-D modeling. G. Ramanujam and A. D. Koussis C
C
,t
~5 4
/
/
|
,,'N
,.¢
L
C,
##
/
I
I I
i /
I
~ 7
?'\
,/,"
//i
~
k.~,~'i
~..-,~ / o
,5
le,
2 J-
C,
~
,,o
o"
/. ,,/i/-~ /
I I
/
01
o
0.5
I
0
I [
,
t
I
7:,.0
0.5
I
Fig. 2. Grid parameter .feasibility diagram for step 1 oJ the FT S scheme
Fig. 3. Grid parameter fi, asiDility dia.qram ,lot step 2 of the F T S scheme
in which p and q are spatial frequencies. Stability is indicated for IGI~
C ~< 2(P -2 + 2) 1/2 - 2P-L. However, equation (26) is also obtained from the discrete perturbation method by requiring' positivity' of solution. Equation (26) is used to delineate the feasibility region shown in Fig. 3. The largest C's permitted for a range of prescribed P values are shown below; in the C versus 0 diagram of Fig. 3 these points plot as a sloping line because, by equation (16), 0(C, P). The condition (20), from the first half step, also ensures positivity of the corrected longitudinal dispersion coefficient used in the second half step.
max
G = [(2 - C - 2K0 + (C + 2p.~) cos(p6x) - \ f ~ - 1C sin(p6x)]/2
(23)
is the equation of an ellipse centered at (2 - C - 21L<)/2on the real axis, with principal axes lengths of (C + 21,L~)and C (Roache 3, p. 44-45). Stability is indicated when the ellipse lies entirely inside the unit circle, i.e., when I G I ~< 1, which requires that I C + 2~ttx[ <~2,
(24a)
C ~< 2.
(24b)
and
Upon substitution of (17b) for ll.~ and rearrangement, equation (24a) yields - 4 < ~ C 2 + 4(C/P) ~< 4
(25)
The left hand inequality is trivially satisfied. The right hand inequality prescribes the stability limit C ~ 2(P -2 + 1) l'z - 2 P - 1
(26)
which is more restrictive than C ~< 2. Application of the more intuitive, but less general, method of stability analysis by discrete perturbation 3 yields the weaker 'dynamic stability' condition
P C
3 1.24
5 1.64
10 m 1.81 2.0
Compared to the LOD scheme, the FTS scheme's area of feasible grid parameters, Fig. 3, is smaller (LOD a r e a = 1/2; FTS a r e a ~ 1/3) thus imposing a more restrictive grid design.
M O D E L T E S T I N G AND C O M P A R I S O N O F RESULTS The FTS and LOD schemes were applied to a set of 1-D and 2-D problems. In this section we present and compare some results of the two tests. A more complete account can be found in Ref. 18, where the graduate research work of D. Syriopoulou 19 and G. Ramanujam z° is summarized. In the I-D test, a Gaussian concentration hill is transported in uniform flow, with u = 0 . 5 , in an unbounded domain with Dx= 10, 20, 50. The initial concentration field is defined as C(x, O) = Co(x) = e x p [ - ( x - Xo)'-'/(2a0)]
(27)
where a o = s t a n d a r d deviation of the concentration field = 264, and x o = center of mass of the concentration field = 2000. In addition, C(x, t) -> 0, as I xl -> ~ .
(28)
This problem was first solved by the LOD method.
Adv. Water Resources, 1991, Vol. 14, No. 4
187
Comparison ~1 two splitting, alyorithms ]or 2-D modeling. G. Ramanujam and ,4. D. Kous'sis -~ S ' C
---5.,~,C'
-
=7E
= '=.
~
~-
e~" t" a :
5"
-'"
4
P 3 x -
i
2
-,,'--"x 2
z¢x
p
=
-
03 b 02
•
4-
t+ /
~
~
*
.
k
;0
\ 0
~is:ance
2
4
( 7r'o-Sand8',
8
6
©istance
10
12
(T~qousanOs)
Fig. 4. Transport of Gaussian hill in unidirectional flow: Comparison of FTS and gOD s c h e m e s . P = 2
Fig. 5. Transport q/'Gaussian hill in unidirectiona/ flow. Comparison o f FTS aqnd LOD schemes. P = 10
When the FTS scheme was applied to the problem, results improved noticeably. All tests were run with 6x = 200, yielding P ranging from 2 to 10. Some of these comparisons are shown for t = 9 6 0 0 in Figs. 4 to 6. The FTS scheme is seen to give more accurate results for all values of the P with comparable C's. The same problem has been run for P = 5 0 , and the results are equally accurate. The grid design restrictions, however, become very tight for large P; and the limiting case P --, co (pure advection) is modeled exactly only for parameter choices that enable tracking of the characteristics: 0 = 1 and C = 2 for the FTS scheme and 0 = 0 . 5 and C = I for the L O D scheme• The second problem considers solute transport in a domain of semi-infinite width and length• The flow is steady and uniform, and the longitudinal and transverse dispersivities are constant, consequently, the dispersion coefficients are constant• The physical parameters for the system are: u = lm/day; D x = 10m-'/day; D~. = 5m-'/day. The uniformly distributed source is located on the entry boundary of the domain• The initial condition is zero concentration over the entire domain• Since the numerical algorithms require a finite width, a fictitious boundary is
placed sufficiently far from the source (130m) for the time frame of simulation, and the concentration on that boundary is specified to be zero (Dirichlet condition). The other transverse boundary passes through the centre of the source and is a symmetry line on which Neumann conditions apply, as shown in Fig. 7. Longitudinal concentration profiles at y = 0 , for various times, were obtained analytically-' ~and compared to results of the FTS scheme ~9 for several grid parameters. The FTS results are in good agreement with the analytical solution. Table 1 gives the computational parameters for cases 4, 9, 12 and 17, which produced the best results. Figure 8 shows the results for case 9. Computations with C = 1.2, not shown, indicate only a slight diminution of accuracy. Plots of nodal error profiles for cases 4 and 9 at 360 and 720 time units are presented in Figs. 9 and 10. The rather large relative errors towards the edge of the plume are not alarming because they are associated with minute concentration values; thus small discrepancies between the analytical and numerical solutions result easily in large relative errors.
CO,qCen[
--
"3liOn
=~x 32,:
C=0 ~- ~TS(C.:5
~ LOE(C-q
,__,40 -4Z ~ - -
P
:
<
5
k
O:
C,
~ ' 2 g
1.0
Table 1.
2
4 [')tS:~Ce
6
8
:0
]2
(Tr,3usaqds)
Fig. 5. Transport of Gaussian hill in unidirectional flow : Comparison of FTS and LOD schemes, P = 5
188
.
a
Fig. 7. Schematic with specified boundary conditions for two-dimensional transport test problem
4-
C
o'C,'~
Adv. Water Resources, 1991, Vol. 14, No. 4
Computational parameters for the FTS scheme
Case
P
C
6t
6x
63'
4 9 12 17
5 5 2.5 2.5
0.6 0.6 0.4 0.4
30 30 10 10
50 50 25 25
10 5 l0 5
Comparison of two splitting algorithms/'or 2-D modeling." G. Ramanujam and A. D. Koussis Concentra, tion
1
7-able 2.
!~
081
---~
A n l I y I I c solution
"
Case 9 (FTS)
Computational parameters liar the FTS scheme
Case
P
C
5 6 7
5 5 5
0.6 0.6 0.4
ot 30 30 20
Jx 50 5O 50
ov 10 5 l0
Concentration ~,~
ol
'"L200° "47 °4 0,"0 :C;
0
600
...... 1000
800
--
An&lytl¢ solution
•
0.8
C I I t 7 (LOO)
1200
Distance 0.6
Fig. 8. Comparison of FTS and analytic solutions for the two-dimensional transport test problenl
0.4
-...
Nodal Error (%) 2O
0.2 - - - Case 4
t'180
- - ~ Case 9
t'360
t'720
t
I 1
15
200
400
soo
8oo
looo
1200
Distance
10
//: t
i ~OO
-5
0
s 200
I 300
i 4 O0
_ _ . a _ _
SOt[,
__
Fig. 11. Comparison of LOD and analytic solutions for the two-dinlensional transport test problenl
~
. . . . . . . . .
t~i"t I
T 0C
Distance
Fig. 9. Nodal errors of FTS scheme solution at 360 days for the two-dimensional transport test problem N o d a l [ rr,.,t t'~,) 40 ................................ -
*
(,~l:.e~ 4
--)
"
Case
g
30
20
720 Days 10
o 0
short times. Figures 12 and 13 show error profiles along the plume centreline ( y = 0) at 360 and 720 time units. The comparison of the error profiles shows that, over the range of test conditions, the FTS scheme has a slight edge in accuracy over the LOD scheme. Since both schemes are formally second order accurate in space and time, this is credited to the use of the most recent results and of a two-dimensional spatial operator in each half-step of the FTS scheme. However, with the specified grid parameters, the integration by the FTS scheme requires about 20% more time on a PC/AT than the LOD scheme. The efficiency of the FTS computation can be improved by using the permitted (and often required) larger C's relative to the LOD scheme for the same P, but at a slight accuracy penalty.
__.~:_~..@... -- -i 200
400
600 Distance
800
1000
~200
Noda: Error (%j
Fig. 10. Nodal errors of FTS scheme solution at 720 days for the two-dimensional transport test problem -lO
Syriopoulou L9 has solved the same problem with the LOD algorithm for different sets of computational parameters and also obtained very good agreement with the analytical solution. Table 2 shows the computational parameters for some of the cases that produced excellent results. Figure 11 shows the L O D results for case 5 at 180, 360, 560 and 720 time units. The agreement with the analytical solution is very good with an increased, though still acceptable, discrepancy at the edge of the plume for
seo
t
-20
i/
- 80
-"
0
CJse 7
lO0
200
900
,aOO
500
600
700
Distance
Fig. 12. Nodal errors of LOD scheme solution at 360 days for the two-dimensional transport test problem
Adv. Water Resources, 1991, Vol. 14. No. 4
189
Comparison oJ" two splitting algorithms f o r 2-D modeling: G. Ramanujam and A. D. Koussis
the strict constraints imposed by equation (26) and increases the competitiveness of the FTS scheme; the degree of relaxation increases with increasing Dy, which in areal models may be one tenth of D~ or greater. The imposed restrictions may be also increasingly relaxed the smoother the input functions are.
x\
-lO
720 Days -20
CONCLUSIONS - 30
Case 7 -40
2GO
400
600
800
1000
1200
Distance Fig. 13. Nodal errors o f L O D scheme sohttion at 720 days f o r the two-dimensional transport test problem
GRID DESIGN CONSIDERATIONS Whether the often marginal gain in accuracy achievable by the computationally more demanding FTS scheme is justified, is an issue that should be decided on a case by case basis. From a practical viewpoint, one important aspect in such a decision is grid design. This aspect can be explored with the help of the feasilibity diagrams. The range of feasible grid parameters is more restricted for the FTS scheme than for the L O D scheme, However, the comparison is not straightforward (at least not in a visual sense) because the match condition in the FTS scheme is 0 = f ( P , C ) while in the L O D it is 0 = f ( P ) . The feasibility regions for both schemes are examined in Table 3 in the range 2~
In this comparative study, care was taken to set the two methods on equal footing: both models compute longitudinal transport in the first step by the Matched Artificial Method scheme and both are second order accurate in space and time. In accord with expectations, the FTS form is more accurate than the simpler L O D procedure. However, in the cases tested the gain is slight and must be paid by an efficiency penalty. When comparable grid parameters are used, the FTS scheme is computationally more expensive than the L O D scheme by about 20%. An advantage of the FTS method is that it allows, and in fact requires, larger C's for accuracy. It is thus possible for the FTS scheme to equal, or even exceed, the efficiency of the L O D scheme. Decisive in the choice between competing numerical methodologies, such as the L O D and FTS/ADI schemes, is often the flexibility the modeller is afforded in grid design. This is especially important in accommodating the demands of field applications. In this regard, the smaller grid parameter range for physically feasible results permitted by the FTS scheme when transverse dispersion is very weak must be viewed as a drawback.
ACKNOWLEDGEMENTS We wish to acknowledge the support received for this study from the United States Geological Survey, Grant # 14-08-0001-G1296, and from Vanderbilt University.
REFERENCES 1 2 3 4 5 6
World Congress International Association of Mathematics and Computers #~ Simulation, Concordia University, Montreal,
Table 3. Comparison offeasible regionsfor the LO D and FTS schemes
P 2 5 10 :c
190
Range of C's LOD 0.0 - 2.0 0.6- 1.4 0.8 - 1.2
7 FTS 0.0 - 1.24 1.2- 1.64 1.6- 1.81
I
Adv. Water Resources, 1991, Vol. 14, No, 4
2
Yanenko,N. N. The MethodofFractional Steps, Springer-Verlag, 1971 Lapidus,L. and Pinder, G. F. Numerical Solution of Partial Differential Equations in Science and Eng#wering, J. Wiley & Sons, New York, 1982 Roache,P. J. Computational Fluid Dynamics, Revised Printing, Hermosa Publishers, A[burqueque, NM, 1976 Cunge,J. A. On the subject of a flood propagation computation method (Muskingum method), J. Hydr. Research, IAHR, 7(2) 1967 205-230 Syriopoulou,D., and Koussis,A. D. Two-dimensionalmodelling of advection-dominated solute transport in groundwater, Hydrosoft, 1988 1(2)63-70 Frind,E. O. and Pinder, G. F. The principal direction technique for solution of the advection-dispersion equation, Proc. Tenth
8 9
Canada, 1982 Frind, E. O., Matanga, G. B. and Cherry, J. A. The dual formulation of flow for contaminant transport modeling. 2: The Borden aquifer, Water Resour. Res., 1985, 21(2) 170-182 Bear, J. Dynamics of Fluids in Porous Media, American Elsevier, 1972 Frind, E. O. and Hokkanen, G. E. Simulation of the Borden Plume using the alternating direction Galerkin technique, Water Resour. Res., 1987 23(5) 918-930
C o m p a r i s o n o f two splitting algorithms J o t 2-D m o d e l i n g It)
11 12
Syriopoulou. D. and Koussis A. D. 2-D Modeling of ad~ection-dominated transport in groundv.ater b? the matched artificial dispersivity method. Water Resour. Res. 1991. 27151 865-872 Hirt. C. W., Heuristic stability theory for finite-difrerence equations. J. Computational Physics, 1968, 2, 339-355 Frind. E. O. The principal direction technique for adjectivedispersive transport simulation in three dimensions. Proc. Fil?h
17
18
13
14 15
16
Bov, en, J. D., Koussis. A. D. and Zimmcr, D. T. Storm drain design: DitTusi~e v~a~e model for PCs. d. Hvdraul. Eml., .4m Soc, Cir. Eml., 1989. 115~8) 1135-1150 Koussis. A. D,. Syriopoulou, D. and Ramanujam. G.. Computation
Of Three-DimeHsiomd .4dvection-Dominated Solute Transport in Saturated Aq,tit~'r~, Report No. USGS 14-08-0001-G1296, NTIS. 19
International Cot{lerence on Finite Elements in Water Revources. Springer-Verlag, New York. 1984. 363-380 Burnett, R. D. and Frind. E. O. Simulation of contaminant transport in three dimensions. 1: The alternating direction Galerkin technique, Water Resour. Res., 1987 23(41 695-705 Peaceman, D. W. Fundartwntals ~1" Numerical Reservoir Simuhttion. Elsevier, 1977 Leonard, B, P., A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comp. Meths. Appl. Mech. Eng., 1979. 19, 59-98 Koussis, A. D., Saenz, M. A. and Tollis, I. G. Pollution routing in streams J. Ht'draul. Eng.. Am. Soc. Cir. Eng., 1983. 109(12) 1636-1651
G. Rarnanujam and .4, D. Koussis
20
1989 Syriopoulou. D,, Two-Dimensional Modellin,," of Cotttanzimtnt Phones in Advccti(m-Dominated Groumtwuter Transport, Ph.D. Disertation. Environmental and Water Resources Engineering. Vanderbilt University, Nashville, Tennessee, 1987 Ramanujam. G..t, hdtidimensiottal Mo&,lin~ of Solute Transport
in Grozmd~aler. Relations of Time Stepping to .qccuracy atut of Dimot.~iomtlitv to Dispersion. M.S. Thesis. Vanderbilt University. 21
1989 Clcary, R. W., Groundwater Pollution and tlydrology, Report 78-WR-15, Water Resources Program, Princeton University. 1978
Ado. W a t e r Resources, 199l, Vol. 14, No. 4
191
Numerical, computation, transport, solute, advection, dispersion, model, groundwater, fractional step, splitting, alternating directions.
INTRODUCTION Fractional steps or splitting methods I were devised to increase efficiency in the numerical solution of multidimensional problems. By treating each spatial dimension implicitly in separate sub-steps, the single, large system of equations of fully implicit methods is decoupled and a series of smaller systems with matrices of smaller bandwidths results that are simpler to solve. In the Peaceman-Rachford Alternating-Direction Implicit (ADI) scheme ~ this is achieved by representing terms in a single dimension implicitly; terms in the remaining dimensions are approximated explicitly. Typically, the solution is advanced over one time step in sub-steps according to the dimensionality of the problem (6t/2 in 2-D, &/3 in 3-D). Yanenko ~ has also proposed a specialized ADI algorithm, which he terms a splitting scheme, that decomposes the governing equation to a number of strictly one-dimensional equations. Lapidus and Pinder-' refer to it as the Local One Dimensional (LOD) approach. We report on the comparison between two related ADI algorithms, one cast in the LOD form and the other in the traditional ADI form in which the space difference operator is multi-dimensional at each fractional step. The LOD scheme is derived in Yanenko 1 and Lapidus and Pinder ~ as a fractional steps method, but it can also be viewed as an approximation in two spatial process-steps, each operating over the entire time interval. Henceforth, we will refer to the standard ADI algorithm as Fractional Time Stepping or FTS method in distinction to the LOD methodology. Both algorithms employ an adaptation of Paper accepted January 1991. Discussion closes February 1992.
© 1991 Elsevier Science Publishers Ltd.
the 'pseudo-viscosity' technique 3, due to Cunge 4, which we term the Matched Artificial Dispersivity (MAD) scheme. Having developed a 2-D LOD groundwater solute transport model 5, we now explore whether retaining both spatial dimensions in each sub-step improves accuracy, and if so, at what cost and under what conditions. The algorithms are designed to integrate the mass balance equation along the principal directions of transport 6. These are the streamlines and a set of corresponding orthogonal coordinates (in media with isotropic hydraulic conductivity, the equipotentials). In this formulation, the transport equation has a single advective term and no cross-derivative terms. There is thus, principally, no difference in computing on Cartesian or on curvi-linear grids. Additionally, transverse numerical diffusion due to non-alignment of grid and flow lines is eliminated. Because the cost of generating the grid nodes on the curvi-linear coordinate system is substantial, the principal directions approach is limited in practice to steady flow fields. It is further assumed that the coordinates' curvature is small. Our interest here is in solute transport in groundwater, where the assumption of steady state flow is common practice 7. G O V E R N I N G EQUATION In many applications, transport of solutes in groundwater can be modelled adequately in two dimensions. With the principal directions of transport as coordinates, the mass balance equation for a conservative dissolved chemical
Adv. Water Resources, 1991, Vol. 14, No. 4
183
Comparison o f two splitting algorithms f o r 2-D modeling. G. Ramanujam and A. D. Koussis
species reads s ?C ?t + u?C.?x = u[?(x,?C/?x)/?x + ((:~.?C.?y)/(y] (1
in which C is the solute concentration, u is the linear pore velocity, x is the curvi-linear coordinate in the direction of flow, y is the coordinate transverse to the flow direction, :z., and :ty are the respective longitudinal and transverse dispersivities, and t is time. If the dispersivities vary so slowly in space (s) that terms (?:t//Ps) (?C/?s) are small compared to Di~ 2 C/~s 2, where D/= ucq, (1) takes the form £C/(t + u?C/?x = Dx 72C/?x2 -}- DyO2C/?y 2
(2)
The longitudinal and transverse dispersion coefficients, D., and D~., are not constant as u = u(x, y), in general. Dispersivities, however, are usually calibrated as constants (as in Ref. 9) because of uncertainty about the field conditions (geometry, flow field, source and observations). Furthermore, to describe the transport of solutes undergoing linear, reversible, equilibrium absorption, u is divided by a 'retardation factor 's. The numerical developments are based on equation (2) with constant coefficients. Space variable u, Dx and D~. can be accommodated by local averaging over grid elements. However, constant grid spacings are implied in all references to truncation errors.
advection dominates, the solution of the pure advection equation is a first approximation to the solution of the complete parabolic equation in the flow direction. On this premise, Cunge a develops an explicit, two-point finite difference (FD) approximation of the advection equation by modelling advection by a time-centred, first order-inspace accurate FD approximation, and by off-centring the time derivative in space (weighting coefficient 0 < 0.5). The discrete model of (3a), with L,. = - ua.,, thus reads ( o / , s t ) ( c * _ ,.~ - c 7 - ~ . j )
+ [(1 - o ) / ~ t ] ( c * j
+ (, 2`5x) [(c,*4 - c ~ _ ~ . j ) + ( c * j - c *
- c,*j) ~4)] = 0
where 3x is the longitudinal increment, ,5t is the time increment, and i, j are space indices, and * indicates the preliminary solution after one ,5t due to longitudinal transport. Rearrangement yields the following directionally explicit FD scheme, in which C = u,st/,sx is the Courant number, C*.=
K , C * _ ~ . j + K , _C ~ _ t . j
+ K3C*
Following the basic concept of the LOD method, equation (2) is split into two one-dimensional equations, here over one time step: (C* - C")/bt = LxC';
(3a)
and
K~ = (C + 20)/(2 + C - 20)
(5a)
Kz = (C - 20)/(2 + C - 20)
(5b) (5c)
The truncation error of this approximation gives rise to numerical diffusion. Cunge 'tunes' this numerical diffusion, to account for the initially neglected longitudinal dispersion term, and obtains second order accuracy in space. He accomplishes this by matching the numerical diffusion coefficient to the dispersion coefficient, the former being derived from a Taylor series expansion to second order. The match specifies the value of the spatial weighting coefficient 0 in terms of the grid Peclet number, P = u6x/D~:
( C n + 1 _ C * ) / 6 t ~- L y C ~
(3b)
where L x = - u6 x + Dx6 z, L~. = Oy62 and ,5i is a first order and ,52 is a second order difference operator. C* is a preliminary solution at the end of one time step ,5t due to longitudinal transport alone. This result is then used as initial condition for the computation of lateral dispersion, again over one 6t. The spatial derivatives are evaluated at some time level n ~< ~ ~
184
(5)
where
K = (2 - C - 20)/(2 + C - 20)
T H E LOD S C H E M E
(4)
Adc. Water Resources, 1991, Vol. 14, No. 4
0 = 0.5 - P - 1
(6)
The scheme is unconditionally stable for physically realizable processes, i.e. those with positive diffusion/ dispersion ~t. The Courant and grid Peclet numbers are subject to mild constraints that derive from the requirement that the scheme's coefficients be positive. Unphysical oscillations are thus eliminated and physically reasonable results ensured'*'5; this control of numerical dispersion does not interfere with the scheme's ability to compute steep fronts. Because of the time-centring of advection, and hence of dispersion (interpreting C* as C "+~ from longitudinal transport, or as C "+~/" in the FTS format), the scheme is also second order accurate in time. In the transverse step, second order accuracy is obtained by centring the second derivative in space and time 1, i.e. (C~.~ ~ - C*j)/6t = Dy(`5~I/.;C "+1 + 6~1,,~C*)/2
(7)
in which 62C ] 14 = (Ci.j+ t - 2C/,j + C~.j_ 1 )//~j2
(8)
Comparison o f two .splitting algorithms j b r 2-D modeling: G. Ramanujam and A. D. Koussis
and 6y is the transverse space increment. As noted, the analysis of truncation error is simplified when starred values are interpreted as concentrations at t = {n + 1/2),3t. THE FRACTIONAL TIME STEPPING ALGORITHM To select an FTS algorithm in which to incorporate the MAD scheme, we considered a number of alternatives presented by Frind 1-' and by Burnett and Frind t3. We decided on the Peaceman-Rachford scheme t't'~ because of its second order accuracy and simple structure. This scheme consists of two steps Step 1
(C,+ 12 _ C")'(5t/2) = LxC "+ 1;_,+ L~.C"
(9)
Step 2(C "+1 - C " + I ~ 2 ) ( 6 t / 2 ) = L x C " + I ! Z + L ~ . C "+~
(10)
We modify this scheme to employ a special form of the MAD scheme in the first half-step in which advection is evaluated at t "+ ~ia. And in further contrast to the LOD method, we retain and evaluate transverse dispersion at t". Evaluation of the space derivatives is thus centred in time over the first half-step (at t "+ t4), which gives second order time accuracy. The temporal derivative is replaced by a weighted sum (weighting coefficient 0) of two FD approximations of the time derivative over &/2 at nodes i - l , j and i,j. The second sub-step is also centred in time, at /,,+3 4. The scheme's equations are Step 1: L_
l/z = lr. ~ , , ~/-r~,i - +1,j
, 1,j) Ci-
+(1 - O ~, (, _C, ,"+ j i'e
-
(il)
cT, j)]/(&12)
Step 2: L ~ C "+ "'- + L y e "+~ = C,")~ - C','.~:
1/2,/,(}
'9
~, t t l - )
(12)
where L.~=-uax:
for step l,
Ly=Dy6~
(13a)
and Lx= -uf.~+Dx6j;
Ly=Dy6~
for step 2. (13b)
Clearly, the solution at each fractional step involves more effort than in the LOD scheme because the difference equations contain terms from derivative approximations in the x- and )'-directions. But for the same reasons the FTS scheme is expected to be more accurate. After collecting terms, we derive for step 1 C~' ~ x,2
I
n
..t
n
t
= K , C , _ t.j + K;C'~+~.') + K3Ci.j -f- KaC'/.j+ ~ + K'sC'~.~-,
(14)
where K't = 2 0 / [ 2 ( 1 - 0) + C ]
(14a)
K2 = (C - 2 0 ) / [ 2 ( 1 - 0) + C ]
(14b)
K ; = [2(1 - 0 ) -
(14c)
and (14e)
& = D~.&;(dy):
Because of the implementation of MAD, the solution of the first step in this formulation is explicit, and thus simple. In deriving equation (14), longitudinal dispersion was neglected. Yet, upon expanding the grid functions in equation (14) in Taylor series around (xi, yj, t "+ ~s-,) to second order, we obtain a transport equation with the numerical longitudinal diffusion coefficient D v = udx/2 - udxO + ua~SU4
To account for longitudinal dispersion, the numerical diffusion coefficient is, again, matched with the physical dispersion coefficient, D., = D s. As before, the match is achieved through the selection of 0; however, in contrast to the LOD match condition (6), in the FTS method 0 depends on P and on C, 0=0.5-
P-~ + 0.25C
K;, = K ; = i , , , / [ 2 ( 1 - 0) + C ]
(14d)
(16)
As in the L O D scheme, only two points are employed implicitly in the longitudinal direction in the first step of the FTS method. Since the second derivative is modelled as an adjusted truncation error of the advection equation, no downstream boundary condition is imposed, nor can one be imposed. This is of little concern in advectiondominated cases, but limits the method's applicability when diffusion is significant. Moreover, it has the advantage that computing time can be economized by limiting the domain longitudinally to the effective transport area at any instant. In the second step, equation (10), the dispersive terms are approximated by the second-order-accurate, spacecentred F D ratios of equation (8). In attempting to resolve the 'smearing' versus overshoot dilemma arising in advection-dominated transport conditions ~4, we experimented with alternative advection schemes. Some were discarded due to deficiencies in accuracy, others because they imposed too restrictive stability constraints (e.g., Leonard's is multi-point, asymmetric Q U I C K scheme). We decided on a 'diffusion-corrected', two-point upwind advection scheme, which we implemented as follows. Upon substituting in (10) the two-point upwind scheme for the advection and space-centred differences for the dispersion terms, the grid functions are expanded about (xi, yj, t "+ k,2) in Taylor series to second order. We obtain again an advection-dispersion equation, except that the effective longitudinal dispersion coefficient is D.~+O.5uSx-O.25u26t. This distortion is rectified by adding to Dx the terms - (ucSx2) + (u2ijt/4). The corrected coefficient of longitudinal dispersion, (Dx).... = D ~ (uax/2) + (u2Dt/4), nullifies the numerical diffusion of the first order accurate approximation for advection, The second-order accurate scheme for the second step is then: ,+1 1 + (2 + 2/ly)C~,;: 1 - ~'~.'-i,j+ ,, r , + l z -l~yCi.~t-,n+
& . l / J 2 ( 1 - 0) + C ]
(15)
1,2
= (2 - C - 2 # . < ) C ' ~ ~"-' + s i . , ~ + 1,j
+ (C . . .l.l x J ( ~
+ lie
i - 1, j
Adv. Water Resources, 1991, Vol. 14, No. 4
(17)
185
Comparison o f two splitting algorithms Jbr 2-D modeling. G. Ramanujam and ,4. D. Koussis
where :tx = (DO ..... 6t/dx 2 = (D x - 0.5udx + 0.25ue 6t)&/6x 2
= ( C P ) - 0.5C + 0.25C 2
(17a)
and as before, ll,. = Oy(t/~) '2, (14e). Because in (1"}) only the transverse dispersion term is implicit, the problem is effectively one-dimensional. The coefficient matrix of the system of linear equations is tri-diagonal and can be solved efficiently by the Thomas algorithm. Accuracy and Stability
In the previous section, both schemes have been shown to be second order accurate in space and time. This has been ensured by time-centring, by lonigitudinal dispersion matching or correction, and by direct second order approximation of the transverse dispersion terms. This section addresses grid restrictions dictated by stability and by physical feasibility of solution. Physically feasible behaviour imposes positivity on the schemes' coefficients, to ensure that positive inputs cause positive responses. This is necessary for the computationally explicit sub-steps; in fact, conditions more restrictive than those imposed by stability may result. The constraints are summarized in diagrams delineating grid parameter regions that ensure meaningful results, which we term feasibility diagrams and feasibility regions, respectively. Both steps of the LOD scheme are unconditionally stable. The stability of the first step is ensured simply by D~. > 0, or 0 < 0.5, which is automatically satisfied for Dx > 0. And the second step has the unconditionally stable Crank-Nicolson format. Only a physical feasibility analysis is therefore required. Restrictions on the feasible region are imposed solely by the first step and from the requirements that K., > 0 and K 3 > 0 in equation (6) (K ~ > 0 always) 4. They are displayed in Fig. 1, along with the stability restriction 0 ~ 0.5 and the limitation imposed by advection-dominated behaviour P/> 2. In the FTS scheme both steps impose limits, Physical feasibility dictates constraints for step 1, numerical
C
C :
2 :~/, ~p_2~i.
0~<0.5+C4
2-20
=o
'
t N~
DN
> o - 05
0
(18)
C~<40
(19)
For equation (14) to give meaningful results, K'~, K ~ . . . K'5 should not be negative. Substitution of equation (16) for 0 in the denominator of these coefficients yields i-1+(2/P)+(C/2)], which is always positive. The numerators are analyzed to determine the sign of K'~, K; • " K'5. K'~ is positive because the minimum value of 0 is C/4 (when P = 2) and C > 0. K~. and K'5 are positive because D~. is positive. To show that K', should not be negative, we consider a zero initial condition and a positive input to the system. Since the downstream concentration response must be positive and only the second term on the right hand side of equation (14) remains, K"2 must be positive, whence c > 20.
(20)
or
0~1
(21)
The conditions (19)--(21) are most restrictive, and thus define the feasibility region for the first half-step shown in Fig. 2. The second step of the FTS method is conditionally stable. When the stability of equation (17) is analyzed by the von Neumann method, the amplification factor 3 obtained is G = [2 - C - 2px + (C + 2p x) cos(p6x)
05
Fi 9. 1. Grid parameter feasibifity diagram f o r the L O D scheme
186
C1>40-2
In fact with the equal sign, equation (18) applies to pure advection, P--+ .m, the correct solution of which is obtained for C = 2, 0 = 1. The scheme is thus stable for any P. That transverse dispersion enhances stability follows from the amplification factor of the yon Neumann stability analysis, which is maximized when Dy is eliminated (by choosing the y-spatial frequency that minimizes the denominator). The feasibility constraints for step 1 are obtained by examining the coefficients in equation (14) and the matching condition (16). Experience with the MAD method ~6"~7 suggests that P = 2 is a meaningful lower limit to characterize advection domination. Application of P > 2 to (16) yields
2(-0)>0, femible
or
For a negative step input from 1 to 0, after the input has ceased only the third, fourth and fifth terms on the right side of equation (14) remain, the last two of which are always positive. On a grid with a sufficiently high transverse resolution we may assume that C~'+~.j~ C7+ ~,j+ ~ ~ C'/+ ~.j_ ~ = 1, which effectively eliminates + 1/9 transverse dispersion. We then have for C ~ l~+ ~.'f > 0
ib [ ,
stability for step 2. In step 1, stability is satisfied automatically through the matching D x = D x, but 0=0(C, P). This can be seen from the positivity of physical dispersion coefficients, D., = D.v > 0, which can be interpreted as a stability constraint~: with equation (15), it yields
Adv. Water Resources, 1991, Vol. 14, No. 4
- \ f ~ - 1C sin(p6x)]/[2 + 2lty(l - cos(q6y)] (22)
Comparison
of two splitting algorithms jor 2-D modeling. G. Ramanujam and A. D. Koussis C
C
,t
~5 4
/
/
|
,,'N
,.¢
L
C,
##
/
I
I I
i /
I
~ 7
?'\
,/,"
//i
~
k.~,~'i
~..-,~ / o
,5
le,
2 J-
C,
~
,,o
o"
/. ,,/i/-~ /
I I
/
01
o
0.5
I
0
I [
,
t
I
7:,.0
0.5
I
Fig. 2. Grid parameter .feasibility diagram for step 1 oJ the FT S scheme
Fig. 3. Grid parameter fi, asiDility dia.qram ,lot step 2 of the F T S scheme
in which p and q are spatial frequencies. Stability is indicated for IGI~
C ~< 2(P -2 + 2) 1/2 - 2P-L. However, equation (26) is also obtained from the discrete perturbation method by requiring' positivity' of solution. Equation (26) is used to delineate the feasibility region shown in Fig. 3. The largest C's permitted for a range of prescribed P values are shown below; in the C versus 0 diagram of Fig. 3 these points plot as a sloping line because, by equation (16), 0(C, P). The condition (20), from the first half step, also ensures positivity of the corrected longitudinal dispersion coefficient used in the second half step.
max
G = [(2 - C - 2K0 + (C + 2p.~) cos(p6x) - \ f ~ - 1C sin(p6x)]/2
(23)
is the equation of an ellipse centered at (2 - C - 21L<)/2on the real axis, with principal axes lengths of (C + 21,L~)and C (Roache 3, p. 44-45). Stability is indicated when the ellipse lies entirely inside the unit circle, i.e., when I G I ~< 1, which requires that I C + 2~ttx[ <~2,
(24a)
C ~< 2.
(24b)
and
Upon substitution of (17b) for ll.~ and rearrangement, equation (24a) yields - 4 < ~ C 2 + 4(C/P) ~< 4
(25)
The left hand inequality is trivially satisfied. The right hand inequality prescribes the stability limit C ~ 2(P -2 + 1) l'z - 2 P - 1
(26)
which is more restrictive than C ~< 2. Application of the more intuitive, but less general, method of stability analysis by discrete perturbation 3 yields the weaker 'dynamic stability' condition
P C
3 1.24
5 1.64
10 m 1.81 2.0
Compared to the LOD scheme, the FTS scheme's area of feasible grid parameters, Fig. 3, is smaller (LOD a r e a = 1/2; FTS a r e a ~ 1/3) thus imposing a more restrictive grid design.
M O D E L T E S T I N G AND C O M P A R I S O N O F RESULTS The FTS and LOD schemes were applied to a set of 1-D and 2-D problems. In this section we present and compare some results of the two tests. A more complete account can be found in Ref. 18, where the graduate research work of D. Syriopoulou 19 and G. Ramanujam z° is summarized. In the I-D test, a Gaussian concentration hill is transported in uniform flow, with u = 0 . 5 , in an unbounded domain with Dx= 10, 20, 50. The initial concentration field is defined as C(x, O) = Co(x) = e x p [ - ( x - Xo)'-'/(2a0)]
(27)
where a o = s t a n d a r d deviation of the concentration field = 264, and x o = center of mass of the concentration field = 2000. In addition, C(x, t) -> 0, as I xl -> ~ .
(28)
This problem was first solved by the LOD method.
Adv. Water Resources, 1991, Vol. 14, No. 4
187
Comparison ~1 two splitting, alyorithms ]or 2-D modeling. G. Ramanujam and ,4. D. Kous'sis -~ S ' C
---5.,~,C'
-
=7E
= '=.
~
~-
e~" t" a :
5"
-'"
4
P 3 x -
i
2
-,,'--"x 2
z¢x
p
=
-
03 b 02
•
4-
t+ /
~
~
*
.
k
;0
\ 0
~is:ance
2
4
( 7r'o-Sand8',
8
6
©istance
10
12
(T~qousanOs)
Fig. 4. Transport of Gaussian hill in unidirectional flow: Comparison of FTS and gOD s c h e m e s . P = 2
Fig. 5. Transport q/'Gaussian hill in unidirectiona/ flow. Comparison o f FTS aqnd LOD schemes. P = 10
When the FTS scheme was applied to the problem, results improved noticeably. All tests were run with 6x = 200, yielding P ranging from 2 to 10. Some of these comparisons are shown for t = 9 6 0 0 in Figs. 4 to 6. The FTS scheme is seen to give more accurate results for all values of the P with comparable C's. The same problem has been run for P = 5 0 , and the results are equally accurate. The grid design restrictions, however, become very tight for large P; and the limiting case P --, co (pure advection) is modeled exactly only for parameter choices that enable tracking of the characteristics: 0 = 1 and C = 2 for the FTS scheme and 0 = 0 . 5 and C = I for the L O D scheme• The second problem considers solute transport in a domain of semi-infinite width and length• The flow is steady and uniform, and the longitudinal and transverse dispersivities are constant, consequently, the dispersion coefficients are constant• The physical parameters for the system are: u = lm/day; D x = 10m-'/day; D~. = 5m-'/day. The uniformly distributed source is located on the entry boundary of the domain• The initial condition is zero concentration over the entire domain• Since the numerical algorithms require a finite width, a fictitious boundary is
placed sufficiently far from the source (130m) for the time frame of simulation, and the concentration on that boundary is specified to be zero (Dirichlet condition). The other transverse boundary passes through the centre of the source and is a symmetry line on which Neumann conditions apply, as shown in Fig. 7. Longitudinal concentration profiles at y = 0 , for various times, were obtained analytically-' ~and compared to results of the FTS scheme ~9 for several grid parameters. The FTS results are in good agreement with the analytical solution. Table 1 gives the computational parameters for cases 4, 9, 12 and 17, which produced the best results. Figure 8 shows the results for case 9. Computations with C = 1.2, not shown, indicate only a slight diminution of accuracy. Plots of nodal error profiles for cases 4 and 9 at 360 and 720 time units are presented in Figs. 9 and 10. The rather large relative errors towards the edge of the plume are not alarming because they are associated with minute concentration values; thus small discrepancies between the analytical and numerical solutions result easily in large relative errors.
CO,qCen[
--
"3liOn
=~x 32,:
C=0 ~- ~TS(C.:5
~ LOE(C-q
,__,40 -4Z ~ - -
P
:
<
5
k
O:
C,
~ ' 2 g
1.0
Table 1.
2
4 [')tS:~Ce
6
8
:0
]2
(Tr,3usaqds)
Fig. 5. Transport of Gaussian hill in unidirectional flow : Comparison of FTS and LOD schemes, P = 5
188
.
a
Fig. 7. Schematic with specified boundary conditions for two-dimensional transport test problem
4-
C
o'C,'~
Adv. Water Resources, 1991, Vol. 14, No. 4
Computational parameters for the FTS scheme
Case
P
C
6t
6x
63'
4 9 12 17
5 5 2.5 2.5
0.6 0.6 0.4 0.4
30 30 10 10
50 50 25 25
10 5 l0 5
Comparison of two splitting algorithms/'or 2-D modeling." G. Ramanujam and A. D. Koussis Concentra, tion
1
7-able 2.
!~
081
---~
A n l I y I I c solution
"
Case 9 (FTS)
Computational parameters liar the FTS scheme
Case
P
C
5 6 7
5 5 5
0.6 0.6 0.4
ot 30 30 20
Jx 50 5O 50
ov 10 5 l0
Concentration ~,~
ol
'"L200° "47 °4 0,"0 :C;
0
600
...... 1000
800
--
An&lytl¢ solution
•
0.8
C I I t 7 (LOO)
1200
Distance 0.6
Fig. 8. Comparison of FTS and analytic solutions for the two-dimensional transport test problenl
0.4
-...
Nodal Error (%) 2O
0.2 - - - Case 4
t'180
- - ~ Case 9
t'360
t'720
t
I 1
15
200
400
soo
8oo
looo
1200
Distance
10
//: t
i ~OO
-5
0
s 200
I 300
i 4 O0
_ _ . a _ _
SOt[,
__
Fig. 11. Comparison of LOD and analytic solutions for the two-dinlensional transport test problenl
~
. . . . . . . . .
t~i"t I
T 0C
Distance
Fig. 9. Nodal errors of FTS scheme solution at 360 days for the two-dimensional transport test problem N o d a l [ rr,.,t t'~,) 40 ................................ -
*
(,~l:.e~ 4
--)
"
Case
g
30
20
720 Days 10
o 0
short times. Figures 12 and 13 show error profiles along the plume centreline ( y = 0) at 360 and 720 time units. The comparison of the error profiles shows that, over the range of test conditions, the FTS scheme has a slight edge in accuracy over the LOD scheme. Since both schemes are formally second order accurate in space and time, this is credited to the use of the most recent results and of a two-dimensional spatial operator in each half-step of the FTS scheme. However, with the specified grid parameters, the integration by the FTS scheme requires about 20% more time on a PC/AT than the LOD scheme. The efficiency of the FTS computation can be improved by using the permitted (and often required) larger C's relative to the LOD scheme for the same P, but at a slight accuracy penalty.
__.~:_~..@... -- -i 200
400
600 Distance
800
1000
~200
Noda: Error (%j
Fig. 10. Nodal errors of FTS scheme solution at 720 days for the two-dimensional transport test problem -lO
Syriopoulou L9 has solved the same problem with the LOD algorithm for different sets of computational parameters and also obtained very good agreement with the analytical solution. Table 2 shows the computational parameters for some of the cases that produced excellent results. Figure 11 shows the L O D results for case 5 at 180, 360, 560 and 720 time units. The agreement with the analytical solution is very good with an increased, though still acceptable, discrepancy at the edge of the plume for
seo
t
-20
i/
- 80
-"
0
CJse 7
lO0
200
900
,aOO
500
600
700
Distance
Fig. 12. Nodal errors of LOD scheme solution at 360 days for the two-dimensional transport test problem
Adv. Water Resources, 1991, Vol. 14. No. 4
189
Comparison oJ" two splitting algorithms f o r 2-D modeling: G. Ramanujam and A. D. Koussis
the strict constraints imposed by equation (26) and increases the competitiveness of the FTS scheme; the degree of relaxation increases with increasing Dy, which in areal models may be one tenth of D~ or greater. The imposed restrictions may be also increasingly relaxed the smoother the input functions are.
x\
-lO
720 Days -20
CONCLUSIONS - 30
Case 7 -40
2GO
400
600
800
1000
1200
Distance Fig. 13. Nodal errors o f L O D scheme sohttion at 720 days f o r the two-dimensional transport test problem
GRID DESIGN CONSIDERATIONS Whether the often marginal gain in accuracy achievable by the computationally more demanding FTS scheme is justified, is an issue that should be decided on a case by case basis. From a practical viewpoint, one important aspect in such a decision is grid design. This aspect can be explored with the help of the feasilibity diagrams. The range of feasible grid parameters is more restricted for the FTS scheme than for the L O D scheme, However, the comparison is not straightforward (at least not in a visual sense) because the match condition in the FTS scheme is 0 = f ( P , C ) while in the L O D it is 0 = f ( P ) . The feasibility regions for both schemes are examined in Table 3 in the range 2~
In this comparative study, care was taken to set the two methods on equal footing: both models compute longitudinal transport in the first step by the Matched Artificial Method scheme and both are second order accurate in space and time. In accord with expectations, the FTS form is more accurate than the simpler L O D procedure. However, in the cases tested the gain is slight and must be paid by an efficiency penalty. When comparable grid parameters are used, the FTS scheme is computationally more expensive than the L O D scheme by about 20%. An advantage of the FTS method is that it allows, and in fact requires, larger C's for accuracy. It is thus possible for the FTS scheme to equal, or even exceed, the efficiency of the L O D scheme. Decisive in the choice between competing numerical methodologies, such as the L O D and FTS/ADI schemes, is often the flexibility the modeller is afforded in grid design. This is especially important in accommodating the demands of field applications. In this regard, the smaller grid parameter range for physically feasible results permitted by the FTS scheme when transverse dispersion is very weak must be viewed as a drawback.
ACKNOWLEDGEMENTS We wish to acknowledge the support received for this study from the United States Geological Survey, Grant # 14-08-0001-G1296, and from Vanderbilt University.
REFERENCES 1 2 3 4 5 6
World Congress International Association of Mathematics and Computers #~ Simulation, Concordia University, Montreal,
Table 3. Comparison offeasible regionsfor the LO D and FTS schemes
P 2 5 10 :c
190
Range of C's LOD 0.0 - 2.0 0.6- 1.4 0.8 - 1.2
7 FTS 0.0 - 1.24 1.2- 1.64 1.6- 1.81
I
Adv. Water Resources, 1991, Vol. 14, No, 4
2
Yanenko,N. N. The MethodofFractional Steps, Springer-Verlag, 1971 Lapidus,L. and Pinder, G. F. Numerical Solution of Partial Differential Equations in Science and Eng#wering, J. Wiley & Sons, New York, 1982 Roache,P. J. Computational Fluid Dynamics, Revised Printing, Hermosa Publishers, A[burqueque, NM, 1976 Cunge,J. A. On the subject of a flood propagation computation method (Muskingum method), J. Hydr. Research, IAHR, 7(2) 1967 205-230 Syriopoulou,D., and Koussis,A. D. Two-dimensionalmodelling of advection-dominated solute transport in groundwater, Hydrosoft, 1988 1(2)63-70 Frind,E. O. and Pinder, G. F. The principal direction technique for solution of the advection-dispersion equation, Proc. Tenth
8 9
Canada, 1982 Frind, E. O., Matanga, G. B. and Cherry, J. A. The dual formulation of flow for contaminant transport modeling. 2: The Borden aquifer, Water Resour. Res., 1985, 21(2) 170-182 Bear, J. Dynamics of Fluids in Porous Media, American Elsevier, 1972 Frind, E. O. and Hokkanen, G. E. Simulation of the Borden Plume using the alternating direction Galerkin technique, Water Resour. Res., 1987 23(5) 918-930
C o m p a r i s o n o f two splitting algorithms J o t 2-D m o d e l i n g It)
11 12
Syriopoulou. D. and Koussis A. D. 2-D Modeling of ad~ection-dominated transport in groundv.ater b? the matched artificial dispersivity method. Water Resour. Res. 1991. 27151 865-872 Hirt. C. W., Heuristic stability theory for finite-difrerence equations. J. Computational Physics, 1968, 2, 339-355 Frind. E. O. The principal direction technique for adjectivedispersive transport simulation in three dimensions. Proc. Fil?h
17
18
13
14 15
16
Bov, en, J. D., Koussis. A. D. and Zimmcr, D. T. Storm drain design: DitTusi~e v~a~e model for PCs. d. Hvdraul. Eml., .4m Soc, Cir. Eml., 1989. 115~8) 1135-1150 Koussis. A. D,. Syriopoulou, D. and Ramanujam. G.. Computation
Of Three-DimeHsiomd .4dvection-Dominated Solute Transport in Saturated Aq,tit~'r~, Report No. USGS 14-08-0001-G1296, NTIS. 19
International Cot{lerence on Finite Elements in Water Revources. Springer-Verlag, New York. 1984. 363-380 Burnett, R. D. and Frind. E. O. Simulation of contaminant transport in three dimensions. 1: The alternating direction Galerkin technique, Water Resour. Res., 1987 23(41 695-705 Peaceman, D. W. Fundartwntals ~1" Numerical Reservoir Simuhttion. Elsevier, 1977 Leonard, B, P., A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comp. Meths. Appl. Mech. Eng., 1979. 19, 59-98 Koussis, A. D., Saenz, M. A. and Tollis, I. G. Pollution routing in streams J. Ht'draul. Eng.. Am. Soc. Cir. Eng., 1983. 109(12) 1636-1651
G. Rarnanujam and .4, D. Koussis
20
1989 Syriopoulou. D,, Two-Dimensional Modellin,," of Cotttanzimtnt Phones in Advccti(m-Dominated Groumtwuter Transport, Ph.D. Disertation. Environmental and Water Resources Engineering. Vanderbilt University, Nashville, Tennessee, 1987 Ramanujam. G..t, hdtidimensiottal Mo&,lin~ of Solute Transport
in Grozmd~aler. Relations of Time Stepping to .qccuracy atut of Dimot.~iomtlitv to Dispersion. M.S. Thesis. Vanderbilt University. 21
1989 Clcary, R. W., Groundwater Pollution and tlydrology, Report 78-WR-15, Water Resources Program, Princeton University. 1978
Ado. W a t e r Resources, 199l, Vol. 14, No. 4
191