- Email: [email protected]
Simulation of convection and diffusion processes by standard finite difference schemes and by influence schemes
COMPUTER METHODS IN APPLIED MECHANICS NORTH-HOLLAND PUBLISHING COMPANY
AND ENGINEERING
3.5 (1982) 153-168
SIMULATION OF CONVECTION AND DIFFUSION PROCESSES BY STANDARD FINITE DIFFERENCE SCHEMES AND BY INFLUENCE SCHEMES G.D.
STUBLEY,*
Department
of Mechanical
G.D.
RAITHBY,
Engineering
Revised
A.B.
STRONG
and Department of Physics, Ont., N2L 3G1, Canada
Received 16 June 1981 manuscript received 18 February
and K.A. University
WOOLNER
of Waterloo,
Waterloo,
1982
It is helpful to consider the error in a finite difference solution as arising from two sources: the profile error component arises because the assumed profile used in deriving the scheme does not match the exact solution; the operator error component is due to the failure of the finite-difference operator to accurately simulate the convection-diffusion process. Profile error is relatively easy to quantify but to understand the performance of any differencing scheme to be used in a problem where convection is important, or dominant, requires a test for operator error. The first contribution of this paper is to propose such a test. Upstream difference, central difference, and influence schemes are subjected to the operator error test. While the influence schemes perform well, it is shown that they may still lead to substantial errors when velocities and diffusion coefficients change sharply over the solution domain. The second contribution of the paper is to show how errors from these sources can be greatly reduced.
1. Introduction Over the past two decades, finite difference methods have been developed to solve complex fluid flow problems involving heat and mass transfer. While these methods have proven adequate for some problems, for many problems they are still inadequate. For example, using the present generation of computers and standard finite difference schemes, it is often prohibitively expensive to reduce false diffusion [l] to an acceptable level in two-dimensional problems involving strong flow recirculation. The accurate solution of the equations of motion for complex three-dimensional flows is often completely beyond reach. A major factor that has stood in the way of developing better methods, i.e., methods which yield greater accuracy for a given grid resolution, has been the lack of a clear connection between the approximations introduced in the derivation of the finite difference equations and the resulting solution errors. Traditional methods of measuring the accuracy of schemes (Taylor series analysis) have not necessarily provided a correct measure of accuracy or a clear understanding of solution error, especially when convection effects dominate. In an earlier paper [2] the authors presented a framework within which the origin of solution errors could be clearly understood. This was based on a proposal to divide the total *Present
address:
Concord
Scientific
004.S7825/82/0000-0000/$02.75
Corporation,
Downsview,
@ 1982 North-Holland
Ont. M3H 2V2, Canada.
154
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
error into two parts, a profile error and an operator error, which led to the development of a Results obtained using the influence new family of schemes called ‘influence’ schemes. schemes showed reduced solution errors. In the present paper the concept of an operator error is further developed and the importance of being aware of this error is reinforced. The first section reviews the subdivision of solution error into its profile and operator error components. This is followed by a section in which a test for operator error is proposed and the test applied to upwind and central difference schemes. The concepts and principles of the influence schemes are then reviewed and these too are subjected to the operator error test. In the final section the poor when diffusivity and velocities vary sharply over the performance of influence schemes, solution domain, is revealed and improvements are proposed which greatly reduce the error under these conditions.
2. Review of operator error For simplitity, the conservation equation diffused and generated, will be considered governing equation is
The solution to (1) with appropriate subdomain (Fig. 1). It is convenient conditions) by .z(@)=
for a single scalar, @, which may be convected, here. For two-dimensional steady conditions the
boundary conditions, is sought over a rectangular to represent this problem (equation plus boundary
(2)
B.
Fig. 1. Solution
domain
illustrating
the coordinate
system
and nomenclature.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
The finite difference
approximation
to the problem
defined
155
by (2) is
L(~P) = bp .
(3)
A lower case 4 is used to denote that the finite difference solution is different from the exact solution, @, and the subscript indicates that discrete 4 values will be obtained at each grid (or P) point. Eq. (3) represents a set of algebraic equations with one equation for each interior and boundary grid point. The equation for an interior point P, Fig. 1, would normally connect & with &, +w, & and & and the equation for a boundary point would connect & (the boundary value) with one or more interior (or I) points. To understand why & and &= are different, it is helpful to follow the steps in the derivation of (3) from (1) [l, 31. Eq. (1) is first integrated over the control volume in Fig. 1 to obtain an exact integral equation. To convert this to an algebraic equation, the values and gradients of @ at the control volume faces and the interior distribution of source terms must be approximated in terms of the discrete values at the grid points. If the correct profiles are chosen in this approximation, the values of & will correspond to the exact solution at the corresponding points (i.e., & = @). Stated alternatively, QP will exactly satisfy (3). Normally the correct profiles are not chosen so that & and aP will differ. In this case the substitution of @ into (3) would leave a residual 6bp, i.e.,
L(@p) = By subtracting obtained:
bp - Sbp
(4) from
.
(3)
(4) the following
equation
for the solution
error,
ep = 4P - QP, is
L(ep) = 6bp. The L-operator is an algebraic analog of 9 and therefore simulates the convection and diffusion processes. Eq. (5) suggests that the error, e P, is generated by the sources 6bp but that the ultimate magnitudes and distribution of the errors are established by the simulated convection and diffusion of ep. Changing the profile approximations affects both the magnitudes of 6bp and the simulated transport. A change in the profile approximation may bring the assumed and exact profiles into closer agreement, resulting in smaller values of 6bp, but the resulting change in the L-operator may sometimes cause larger solution errors [2]. This provides an explanation for the possibly perplexing conclusion that a reduction in truncation error does not necessarily reduce the solution error [4, 51. In recognition of the separate roles played by the residuals and the operator in determining the solution error, it is useful to split the error into two portions [2], the profile error, epr, and in (2)) were perfectly simulated by the operator error, cop. If the transport (i.e., the Z-operator the L-operator in (5) the Sbp sources would give rise to the profile error, epr, obtained by solving 9(eP’) = Because
66.
(6) is continuous,
the source
function
S6 must be the continuous
representation
of the
156
schemes for convection and dijjiuion
G.D. Stubley et al., Finite difference and influence
discrete residuals, SbP. The transformation from discrete to continuous is not unique but different choices of distribution have been shown to give nearly identical values of ePr [6]; it is simplest to take 6b as uniform over each control volume with the same total source as 6bP. The remaining error, the operator error, cop, is due to the inadequacy of L as a model of 9. Since ep = egp + e$, (5) and (6) yield the following (symbolic) equation for eS, egp = Lml(Sbp) -
[.2?‘(13b)]p.
0)
Higher order finite difference schemes reduce the profile error component of the total error, but how is the operator error, or the simulation of the convection and diffusion processes, affected? To answer this question it is desirable to develop a test procedure for evaluating the operator error for any given finite difference scheme. Such a test is now described; applications of the test should also clarify the concept of operator error.
3. A test for operator error In (5) (and (6)) the values of 6b, (and 86) are generally non-zero throughout the solution domain. However, because both equations are linear, the same solutions could be obtained by superposition of N simpler solutions, where N is the number of grid points; each of the simpler solutions is obtained by setting SbP (and 86) to zero at all grid points except one. To understand the performance of a given finite-difference scheme it is only necessary to solve one such subproblem. Eqs. (5) and (6) are now respectively solved with the one non-zero residual term arbitrarily set to unity. The solution domain for the proposed test problem is shown in Fig. 2. The residual source term, 6b,, is zero everywhere except over the center control volume. The continuous residual source term, 66, is constant over this control volume, as already described. The diffusivity is
Y’t
+ R
0
1
+I
+
+ +
-A L +
+
+
+++
Cl4
R
d=I
2R
Fig. 2. Layout for the operator test problem.
x+
157
G.D. Stubley et al., Finite difference and influence schemes for convection and dijjiuion
constant and the velocity, also constant, is at an angle 8 to the x* axis. Under (1) which defines the Z-operator in (6) becomes .Z(eP’)=cost3$eP’+sin 6b=
1 0
these conditions
* pr a2epr = s6 2 +-++T
O$epr-$e
ay
@I
if R-$s(x*,y*)sR++, otherwise,
where R defines the solution domain and x” = x/A, y* = y/A, Pe = pVA/I’. The normalized distances and Peclet number are based on the grid-spacing, A, which is constant over the solution domain. On boundaries where there is a convective flux into the solution domain, or for any boundaries where there is a negligible convective flux, the value of ePr is set to zero. For boundaries where there is a convective flux out of the solution domain, the normal gradient of ePr is set to zero. For example, when the velocity is parallel to the x*-axis, the appropriate boundary conditions are epr(O, y *> = 0 ) whereas
for velocity
ePr(O,y *)
epr(x*,O)=O,
angles greater
=0,
$1
=O,
ePr(x*, 2R) = 0 ,
x*=2R
than 0” but less than 90”, the boundary
epr(X*,O)=O,
31
=o, x’=2R
%I aY
conditions
=o. y*=ZR
are
(10)
These boundary conditions have been chosen so that numerical solutions obtained using central difference schemes do not have ‘wiggles’ caused by the outflow boundary conditions as discussed by Roache [3]. The exact solution of (8) subject to these boundary conditions is given in Appendix A. The finite difference approximation of these equations, i.e., of (8) and (9) or (lo), has already been symbolically written as (5) and the solution to the finite difference equation set yields the total solution errors, ep. Since the operator error is simply the additional error caused by this finite difference approximation of 2 by L, e”pP = ep - ePr. Numerical solutions were obtained using two finite difference schemes. A linear profile in 4 between grid points for the convective and diffusive fluxes was assumed, leading to the classical central difference scheme (CDS) for one scheme. A constant stepwise profile was assumed for the convective fluxes leading to the classical upwind difference scheme (UDS) [3] for the other scheme. The exact and numerical results are presented for Peclet numbers (Pe) of 1 and 10 and for flow angles (0) of 0” and 45”. Contours of constant ePr arising from the single source have been evaluated from the exact solution and are plotted as solid lines in Figs. 3-6. Contours of constant e have been evaluated from the finite difference solutions and are plotted as broken lines. The closeness of the broken line to the solid line provides a direct indication of operator error. Since the results in each case are symmetrical about the streamline through the centre of
G.D. Stubley et al., Finite difference and influence schemes for convection and di@siun
1.58
20.07
20.0 Pe-
Ps = 1.0
1.0
8 =45*
8 = 0’
e p’=o.05
15.0 -
i
/
,
/
Y* 10.0-
5.0 -
o.oL 5.0
0.0
IO.0
15.0
0.0
20.0
5.0
I
I
I
10.0
15.0
20.0
X*
X* Fig. 3. CDS and UDS operator test for Pe = 1.0 and t? = 0”.
Fig. 4. CDS and UDS operator test for Pe = 1.0 and 0 = 45”.
the source volume only half of the contours from the numerical solutions are given. The upper halves are the contours based on CDS results and the lower halves are based on UDS results. The contours are for lines of e = 0.05. The results presented in Figs. 3 and 4 are for Pe = 1.0. For this low grid Peclet number the difference between the profile and total error (i.e., operator error) is small for both schemes. The UDS contour is slightly wider than the exact contour for 8 = 4.5”, which indicates that the influence of the residual source is being incorrectly propagated normal to the velocity vector. 20.0
i-
I
Pa = 10.0 8
= 00 IS.0
Y* 10.0
:~~,,,” “.y o.oliI 0.0
5.0
10.0
15.0
20.0
X0
Fig. 5. CDS and UDS operator test for Pe = 10.0 and i3 = 0”.
0.0
5.0
10.0
,
i
15.0
20.0
X+
Fig. 6. CDS and UDS operator test for Pe = 10.0 and 6 = 45”.
G.D. Stubley et al., Finite difference and influence schemes for convection and difision
1.59
The results are much more dramatic for the higher Peclet number, Pe = 10.0. For 0 = O”, the UDS results, shown in Fig. 5, closely match the exact contour except near the source volume. A portion of the CDS results match the exact contour downstream of the source volume but additional closed contours appear upstream of the source. This indicates that the CDS operator incorrectly propagates the influence of the residual source upstream. For 8 = 45” the upstream propagation by the CDS operator is again evident (Fig. 6); negative values of e as large as 63.5% of the value at the centre of the source volume appear. The UDS operator gives a contour that is much wider than the exact contour. This incorrect propagation normal to the velocity vector has been referred to as ‘false’ or ‘numerical diffusion’ [3, 41. This test problem permits the convection-diffusion characteristics of the CDS and UDS operators to be clearly demonstrated and understood. Other numerical operators can also be tested and compared using the same problem.
4. Performance
of influence schemes
When the finite difference equations are derived by conventional means, the choice of the profile shape fixes both the profile error and the operator error. In an attempt to assert a certain amount of separate control over each of the error components, the authors have proposed a new class of finite difference schemes called ‘influence schemes’ [2]. A different motivation led Chen and co-workers [7] to independently propose a similar class of schemes. A brief description of influence schemes and the results of subjecting two such schemes to the test problem described in the previous section will now be given. The P-point in Fig. 7 has, as its immediate neighbours, the 8 points N, NE, E, . . . , NW. An accurate simulation of the convection-diffusion process is desired within the subdomain whose boundary lies through these points (solid lines in Fig. 7) so that the +-value at each boundary point accurately ‘influences’ the value of &. In addition, the presence of sources of 4 within each of the 9 subsections (denoted by dotted lines in Fig. 7) should also accurately influence
l
BOUNDARY SUBDOMAIN /
Fig. 7. Subdomain
for determining
r-
OF
SUBDOMAIN
the influence
coefficients.
160
G.D.
Stubley et al., Finite difference and influence schemes for convection
and di#iuion
the value of &. To obtain accurate influence, the finite difference operator, L, should closely simulate the differential operator in (1). The L-operator was therefore based on an analytical solution to the equation
w [pulp ax
+
“*4+(-j [PVIP-w = rP 2a*++r dX
Pdy2
dY
(11)
r
where the coefficients are constant, evaluated at P, and d1 is constant, but perhaps different, over each of the 9 subsections in Fig. 7. To solve (11) continuous boundary conditions around the subdomain in Fig. 7 are required. This requires a profile assumption to specify a distribution of 4 between the grid points around the subdomain. Linear and quadratic profiles give rise respectively to the linear influence scheme (LIS) and the quadratic influence scheme (01s). The resulting finite difference equation for the P-subdomain has the form
where the & are the values of 4 at the 8 neighbour points, and the the 9 subsections. The analytical solution of (11) subject to the ditions, results in expressions for the ‘influence coefficients’, C1, and coefficients’, DI. The derivation of these expressions, by the Green’s underlines the influence point of view [6]. The expressions for Appendix B. When the two influence schemes, LIS and QIS, are applied to the the last section, with Pe = 1.0 and 13= 0” and 45”, the numerical and
0, are the source terms in prescribed boundary conthe ‘source term influence function method, further C, and Dr are given in test problem, described in analytical solutions (i.e., e
20.0 -
20.0-
Pe = 10.0
PO = 10.0
15.0 -
a. I.S.--eP’=0.05
Y* 10.0-
0.0 0.0
I 5.0
I 10.0
I 15.0
I 20.0
X*
Fig. 8. QIS and LIS operator test for Pe = 10.0 and 6 = 0”.
0.0 / 0.0
I 5.0
I 10.0
1 15.0
I 20.0
X*
Fig. 9. QIS and LIS operator test for Pe = 10.0 and 0 = 45”.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
161
and epr) were virtually identical. For Pe = 10.0, Figs. 8 and 9 show that the QIS leads to for 8 = 0” and 45”, respectively. These figures also show that the LIS accurate predictions operator gives a contour slightly wider than the exact contour for 8 = 45”. It should also be noted that small negative values (3.4% of the value at the centre of the source volume) appear in the QIS solution for Pe = 10.0 and 8 = 45”. These results clearly indicate that the operator errors of the influence schemes have been kept to a low level. While the operator error is principally determined by the operator in (11) the results indicate that the profile approximations still play a slight role in determining the operator error.
5. Improvements
to the influence schemes
Two approximations are required to derive the finite difference equations for the influence schemes, a profile approximation around the subdomain (Fig. 7) and the form of the operator used to obtain the influence coefficients, (11). In the previous section and earlier work [2], a comparison of results based on linear and quadratic profiles has shown that improving the profile approximation reduces solution error. In this section it will be shown that the finite difference operator derived from (11) is still inadequate for some problems, and that a better operator can easily be derived. The finite difference operator for the influence schemes, and Chen’s scheme, is based on an analytical solution of (11) as described. The constant coefficients, assumed in (11) can lead to large errors when the velocity and/or diffusion coefficients vary sharply over the subdomain (Fig. 7). This error can be substantially reduced by modifying the source terms, as will now be shown. Eq. (1) which defines the exact operator, can be written without approximation as
[pfi]p
g
+
[pv’],
a*@ +r”*@ E = r---y 7+d ax dY aY
where
o= d-Cori-[pa],)~-@6-[p6]~)~.
(W
The mass conservation equation has been used to modify the convection terms, the diffusion terms have been expanded, and new apparent mass fluxes have been defined as pti = PU - arlax,
pc = PV - arjay
.
(14)
Eq. (13a) is of the same form as (11) if F is approximated by r, and if 0 is taken as constant over each of the 9 subsections in Fig. 7. Thus expressions for C, and D1 in (12) based on (13a) are identical to those based on (11) except that [piilP and [p& replace [pulp and [pv)lP, respectively. In evaluating d for each subsection, the values of pfi, pv’, a@/dx and a@/ay are evaluated at the grid point associated with that subsection; &?/ax and a@/ay are based on the best available (usually from the last iteration) values of 4 and the derivatives are ap-
162
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
proximated by passing quadratic profiles through the surrounding grid points. The linear and quadratic influence schemes that have been extended in this manner have been named LISE and QISE respectively. The importance of these extensions can be seen by applying the influence schemes to obtain numerical solutions to the following problem, proposed by Baliga and Patankar [8]. The distribution of a conserved scalar @ in a flow between a rotating and fixed cylinder, Fig. 10, is sought on a Cartesian grid within the test section shown in the north-east quadrant. In Cartesian coordinates the problem appears to be two-dimensional with velocity and diffusivity fields which may vary spatially. Analytical expressions for the distributions of the velocity and the scalar @ can be readily obtained by formulating the problem in cylindrical coordinates. The conditions under which solutions were obtained in the present work are: (i) the outer cylinder is fixed and the inner cylinder rotates at an angular velocity wi; (ii) all variables are nondimensionalized by their values at the inner cylinder radius ri, i.e.,
(iii) the diffusivity of the scalar @ is given by r* = 1.0/r*2; (iv) @ is 1.0 at the inner cylinder and 0.0 at the outer cylinder. In cylinder coordinates the expression for the velocity is V,*(r*) =
(r*2 - ro*‘) r*(l.O - rZ2)
and that for the scalar @ is
&.‘I,; Fig. 10. Layout
for rotating
the cylinder
test problem.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
163
(17) For this problem
(1) simplifies
to
where
U
*=
“*yT f3 f%
3
V
*-&
pe=
”
(19)
pirToi’
Numerical solutions to (18) were obtained using the quadratic influence scheme, QIS, the extended influence scheme, QISE, and the exponential difference scheme, EDS [2], which reduces to CDS at low Peclet numbers and to UDS at high Peclet numbers. The exact solution, (17) was used to specify the values on the boundary of the test section. The outer cylinder radius, r i, was set at 3.0 and each scheme was tested for a wide range of Peclet number. The set of finite difference equations, for each test, was solved using the line-by-line method [l], and the solution was considered to have converged when the maximum change of 4 at every grid point was less than 0.001% between iterations. The performance of each scheme was assessed by studying the percentage error in @, that was predicted at the centre of the test section, as the Peclet number was varied. The results, Fig. 11, demonstrate that the solution based on the influence scheme, QIS (without modification), had prohibitively large errors for this problem at small Pe (when diffusive transport dominates). The extended influence scheme, QISE, outperforms, in terms of
VI0
=
0.0
r = l/r’*
0.01
0.1
IO
1.0
100
1000
Pe+ Fig. 11. Midpoint
errors
for the rotating
cylinder
test problem.
164
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
accuracy, both EDS and QIS for all Peclet numbers tested. Similar results were obtained for this problem but under different conditions i.e., where both cylinders rotated at the same angular velocity, or the diffusivity was constant, etc. [6]. These results show the large solutions errors that can arise by using QIS in cases where the velocity and diffusivity vary spatially. By using the QISE, the solution errors can be greatly reduced with little additional numerical effort.
6. Summary This paper makes two separate contributions and, because the basic ideas will be new to most potential readers, a review of some previous work was also presented. First the equations for the error in a finite difference solution were written and the motivation for subdividing the error into ‘profile’ and ‘operator’ components was reviewed. The truncation error of the finite difference scheme provides a qualitative measure of the profile error but no good test for operator error was available in the literature. The first main contribution of this paper was to propose such a test and to subject a central difference scheme, an upwind difference scheme, and two influence schemes to this test; the ideas behind the formulation of influence schemes were reviewed. While influence schemes have desirable properties, it is shown that large solution errors can still exist if the diffusion coefficients and velocities change rapidly over the solution domain. The second contribution of this paper is to show a simple procedure for greatly reducing the solution error for such problems.
Acknowledgment The financial support of the Natural Sciences and Engineering Research Council of Canada, through operating grants and through a scholarship to G.D. Stubley, is gratefully acknowledged.
Appendix A In the text a test for the operator error was discussed. The profile error, solving (8) subject to the boundary conditions given in (9) or (lo), depending the velocity vector. The solution to (8) given in this appendix is obtained function method [6]. The exact solution to (8) subject to the boundary conditions given in (9)
epyx*,y*) = 2
sir@y *) exp[Pe,x
“=l y R + $$ where
*/2 + Pe,y */2] I&
sin(2@?))(cosh(2yR)
+2
sinh(2yR))
epr, is obtained by upon the angle of using the Green’s is
(A4
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
p = nn/2R,
y2 = Pe2/4 + p” ,
Pe,=Pecosf3,
Pe, = Pe sin
165
8,
Rttl2 II
=
I
R-112
evPe,d21 sin(h) dv , sinh(yx*)[cosh(y(2R
- 5)) + %
sinh(y(2R
- &))I ; 6 > **
- x*)) + %
sinh(y(2R
- x*))]
expPe,5/21
d5.
sinh(y5)[cosh(y(2R
;e < x*
I
The exact solution to (8) subject to the boundary conditions given in (10) is identical given in (A.l) except for the eigenvalues which are now given as p = (2n - 1)7r/4R.
Appendix
to that
(A.21
B
The influence coefficients, been obtained [6] by solving expressions are given in this The influence coefficients r = AxlAy
C,..,, etc., and the source term influence coefficients, D,, etc., have (11) analytically using the Green’s function method. The resulting appendix. are functions of three nondimensional parameters:
(grid aspect
ratio),
Pe, = [PU],AX/T~
(local x Peclet
number)
Pe, = [p~]~Ay/r~
(local y Peclet
number).
The resulting
expressions
for each influence -
ci =exdipex+ iPeylmsl 2
coefficient
(-l)“Af,, cosh(a r) m
(B-1)
,
have the form
03.2)
where u’, = a2 + pz, &, = i(2m - 1)~ and a2 = Pez/4r2 + Pet/4. The expressions for the series coefficients, A”, and A”,“, are given in Table B. 1 for LIS and QIS. It should be noted that influence of each corner point is broken into two components so that
G’ = Gvw +
csws
(B.3)
where the expressions for C,, and C,,, have the form given in (B.2). The expressions for the is other influence coefficients are all related to C, and C,,,. A summary of these relationships given in Table B.2. The expressions for the source term influence coefficients have the form
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
166 Table B.1 Expressions
for the series coefficients
Influence coefficient
Boundary condition function
CW
f = {(2’;:,,
o
Am for the influence
Scheme 1
schemes
Series coefficient
LIS where
y = :Pe: + /3$
QIS
f=Y*(2-y*),
CW
Ah
O
C SWW
1
C SW
h =(l-y*)(2-y*)/2,
LIS QIS
k-3+_
Y
peP+&($peLP2) Y Y3
m
o
03.4) where the series coefficient pairs, I&, for each influence coefficient are
The series coefficients, Ii and lj, are given in Table B.3. To aid in the convergence of the series (B.4) for certain coefficients the following identities have been derived.
&
D,
r 0, AY2
=epey"S2 - e Pei4(iPeyS0 + S,) + exp[iPe,
+ tPe,] .z, 2J-~~~$$ m
), m
= eP’~‘4($Pe,(1 - e-Pey’2)So + (1 + e-pey’z)S1)
I, &
D, = e-Pey’2S2+ epPe~4(~Pe,,So- S,) + exp[$Pe, + $Pe,] 2 (- ‘)“lrxr, m=, 2~~ cos(ra,)
03.6)
’
167
G.D. Stubley et al., Finite difference and influence schemes for conuection and diffusion
where 2 SO= Pe$ sinh(Pe,)
coth(Pe,) + &-)(cosh(sPe,)
- cosh(iPe,))
Y
- %sinh(jPe,) + j sinh(iPe,) ,
s, = Pe
siih(Pe Y
) coth(Pe,)(sinh(!Pe,)
+ sinh($r
Y
- i sinh($Pe,) - a cosh(iPe,) , 1 sinh(PeJ2) ” = 2Pe, cosh’(PeJ2)
(3.7)
’
Table B.2 Summary of influence coefficients in terms of functions for C, and C,, C, = Csw(Pex, Pe,, r) + CswPe,, Pex, l/r) -Pe,, l/r) + e-“xC,(PeX. Pe,, r) C,, = Cm(Pex, -Pe,, r) + ewpeyG.&k,, Pe,, I/r) C,, = e-P”YCsww(Pey,-Pe,, l/r) + e-PeW~w(Per, -Pe,,
C, = C,(Pe,, Pe,, r)
c,= Csw(PeY,
Table B.3 Expressions
C, = e-P’xCw(Pe,, Pe,, r) C, = C,(Pe,, -Pe,, l/r)
r)
C, = ewPeyCw(Pey,-Pe,,
for source term influence coefficients
= 2u,,, cosh(ra,)
I PX
e-~ex~2
+
I n
Y
e -3Pe&T
I, = - r/ - c k
sinh $rcr,,,- ra, co& k,,,]
[$e, [$e,
e-3PQl‘i = ry
[-$e,
sinh dru,,, + ru,,, cash sinh
irm,,,
+
?ra,]
ra, cash ircr,,,] - a, Y
&+
4v
=
COS(&~~)[&
F ewPey4 CoS(b> Y
+ iPe,(-
[1p (_I)“(1 2 ey
_
I)“1
e-Pey/Z)
+ pm(l
+ e-*W)]
l/r)
168
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
References [l] S.V. Patankar, Numerical Heat Transfer and Fluid Flow (McGraw-Hill, New York, 1980). [2] G.D. Stubley, G.D. Raithby and A.B. Strong, Proposal for a new discrete method based on an assessment of discretization errors, Numer. Heat Transfer 3 (1980) 411-428. [3] P.J. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, NM, 1972). [4] G.D. Raithby, A critical evaluation of upstream differencing applied to problems involving fluid flow, Comput. Meths. Appl. Mech. Engrg. 9 (1976) 7.5-103. [5] K. Torrance et al., Cavity flows driven by buoyancy and shear, Fluid Mech. 51 (1972) 221-231. [6] G.D. Stubley, A new discrete method for convection-dominated flows based on a clear understanding of solution errors, Ph.D. Thesis, University of Waterloo, Waterloo, Ont., Canada, 1981. [7] C.J. Chen and P. Li, Finite differential method in heat conduction-application of analytical solution technique, ASME Paper 79-WA/HT-50, 1979. [8] B.R. Baliga and S.V. Patankar, A new finite-element formulation for convection-diffusion problems, Numer. Heat Transfer 3 (1980) 395410.
AND ENGINEERING
3.5 (1982) 153-168
SIMULATION OF CONVECTION AND DIFFUSION PROCESSES BY STANDARD FINITE DIFFERENCE SCHEMES AND BY INFLUENCE SCHEMES G.D.
STUBLEY,*
Department
of Mechanical
G.D.
RAITHBY,
Engineering
Revised
A.B.
STRONG
and Department of Physics, Ont., N2L 3G1, Canada
Received 16 June 1981 manuscript received 18 February
and K.A. University
WOOLNER
of Waterloo,
Waterloo,
1982
It is helpful to consider the error in a finite difference solution as arising from two sources: the profile error component arises because the assumed profile used in deriving the scheme does not match the exact solution; the operator error component is due to the failure of the finite-difference operator to accurately simulate the convection-diffusion process. Profile error is relatively easy to quantify but to understand the performance of any differencing scheme to be used in a problem where convection is important, or dominant, requires a test for operator error. The first contribution of this paper is to propose such a test. Upstream difference, central difference, and influence schemes are subjected to the operator error test. While the influence schemes perform well, it is shown that they may still lead to substantial errors when velocities and diffusion coefficients change sharply over the solution domain. The second contribution of the paper is to show how errors from these sources can be greatly reduced.
1. Introduction Over the past two decades, finite difference methods have been developed to solve complex fluid flow problems involving heat and mass transfer. While these methods have proven adequate for some problems, for many problems they are still inadequate. For example, using the present generation of computers and standard finite difference schemes, it is often prohibitively expensive to reduce false diffusion [l] to an acceptable level in two-dimensional problems involving strong flow recirculation. The accurate solution of the equations of motion for complex three-dimensional flows is often completely beyond reach. A major factor that has stood in the way of developing better methods, i.e., methods which yield greater accuracy for a given grid resolution, has been the lack of a clear connection between the approximations introduced in the derivation of the finite difference equations and the resulting solution errors. Traditional methods of measuring the accuracy of schemes (Taylor series analysis) have not necessarily provided a correct measure of accuracy or a clear understanding of solution error, especially when convection effects dominate. In an earlier paper [2] the authors presented a framework within which the origin of solution errors could be clearly understood. This was based on a proposal to divide the total *Present
address:
Concord
Scientific
004.S7825/82/0000-0000/$02.75
Corporation,
Downsview,
@ 1982 North-Holland
Ont. M3H 2V2, Canada.
154
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
error into two parts, a profile error and an operator error, which led to the development of a Results obtained using the influence new family of schemes called ‘influence’ schemes. schemes showed reduced solution errors. In the present paper the concept of an operator error is further developed and the importance of being aware of this error is reinforced. The first section reviews the subdivision of solution error into its profile and operator error components. This is followed by a section in which a test for operator error is proposed and the test applied to upwind and central difference schemes. The concepts and principles of the influence schemes are then reviewed and these too are subjected to the operator error test. In the final section the poor when diffusivity and velocities vary sharply over the performance of influence schemes, solution domain, is revealed and improvements are proposed which greatly reduce the error under these conditions.
2. Review of operator error For simplitity, the conservation equation diffused and generated, will be considered governing equation is
The solution to (1) with appropriate subdomain (Fig. 1). It is convenient conditions) by .z(@)=
for a single scalar, @, which may be convected, here. For two-dimensional steady conditions the
boundary conditions, is sought over a rectangular to represent this problem (equation plus boundary
(2)
B.
Fig. 1. Solution
domain
illustrating
the coordinate
system
and nomenclature.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
The finite difference
approximation
to the problem
defined
155
by (2) is
L(~P) = bp .
(3)
A lower case 4 is used to denote that the finite difference solution is different from the exact solution, @, and the subscript indicates that discrete 4 values will be obtained at each grid (or P) point. Eq. (3) represents a set of algebraic equations with one equation for each interior and boundary grid point. The equation for an interior point P, Fig. 1, would normally connect & with &, +w, & and & and the equation for a boundary point would connect & (the boundary value) with one or more interior (or I) points. To understand why & and &= are different, it is helpful to follow the steps in the derivation of (3) from (1) [l, 31. Eq. (1) is first integrated over the control volume in Fig. 1 to obtain an exact integral equation. To convert this to an algebraic equation, the values and gradients of @ at the control volume faces and the interior distribution of source terms must be approximated in terms of the discrete values at the grid points. If the correct profiles are chosen in this approximation, the values of & will correspond to the exact solution at the corresponding points (i.e., & = @). Stated alternatively, QP will exactly satisfy (3). Normally the correct profiles are not chosen so that & and aP will differ. In this case the substitution of @ into (3) would leave a residual 6bp, i.e.,
L(@p) = By subtracting obtained:
bp - Sbp
(4) from
.
(3)
(4) the following
equation
for the solution
error,
ep = 4P - QP, is
L(ep) = 6bp. The L-operator is an algebraic analog of 9 and therefore simulates the convection and diffusion processes. Eq. (5) suggests that the error, e P, is generated by the sources 6bp but that the ultimate magnitudes and distribution of the errors are established by the simulated convection and diffusion of ep. Changing the profile approximations affects both the magnitudes of 6bp and the simulated transport. A change in the profile approximation may bring the assumed and exact profiles into closer agreement, resulting in smaller values of 6bp, but the resulting change in the L-operator may sometimes cause larger solution errors [2]. This provides an explanation for the possibly perplexing conclusion that a reduction in truncation error does not necessarily reduce the solution error [4, 51. In recognition of the separate roles played by the residuals and the operator in determining the solution error, it is useful to split the error into two portions [2], the profile error, epr, and in (2)) were perfectly simulated by the operator error, cop. If the transport (i.e., the Z-operator the L-operator in (5) the Sbp sources would give rise to the profile error, epr, obtained by solving 9(eP’) = Because
66.
(6) is continuous,
the source
function
S6 must be the continuous
representation
of the
156
schemes for convection and dijjiuion
G.D. Stubley et al., Finite difference and influence
discrete residuals, SbP. The transformation from discrete to continuous is not unique but different choices of distribution have been shown to give nearly identical values of ePr [6]; it is simplest to take 6b as uniform over each control volume with the same total source as 6bP. The remaining error, the operator error, cop, is due to the inadequacy of L as a model of 9. Since ep = egp + e$, (5) and (6) yield the following (symbolic) equation for eS, egp = Lml(Sbp) -
[.2?‘(13b)]p.
0)
Higher order finite difference schemes reduce the profile error component of the total error, but how is the operator error, or the simulation of the convection and diffusion processes, affected? To answer this question it is desirable to develop a test procedure for evaluating the operator error for any given finite difference scheme. Such a test is now described; applications of the test should also clarify the concept of operator error.
3. A test for operator error In (5) (and (6)) the values of 6b, (and 86) are generally non-zero throughout the solution domain. However, because both equations are linear, the same solutions could be obtained by superposition of N simpler solutions, where N is the number of grid points; each of the simpler solutions is obtained by setting SbP (and 86) to zero at all grid points except one. To understand the performance of a given finite-difference scheme it is only necessary to solve one such subproblem. Eqs. (5) and (6) are now respectively solved with the one non-zero residual term arbitrarily set to unity. The solution domain for the proposed test problem is shown in Fig. 2. The residual source term, 6b,, is zero everywhere except over the center control volume. The continuous residual source term, 66, is constant over this control volume, as already described. The diffusivity is
Y’t
+ R
0
1
+I
+
+ +
-A L +
+
+
+++
Cl4
R
d=I
2R
Fig. 2. Layout for the operator test problem.
x+
157
G.D. Stubley et al., Finite difference and influence schemes for convection and dijjiuion
constant and the velocity, also constant, is at an angle 8 to the x* axis. Under (1) which defines the Z-operator in (6) becomes .Z(eP’)=cost3$eP’+sin 6b=
1 0
these conditions
* pr a2epr = s6 2 +-++T
O$epr-$e
ay
@I
if R-$s(x*,y*)sR++, otherwise,
where R defines the solution domain and x” = x/A, y* = y/A, Pe = pVA/I’. The normalized distances and Peclet number are based on the grid-spacing, A, which is constant over the solution domain. On boundaries where there is a convective flux into the solution domain, or for any boundaries where there is a negligible convective flux, the value of ePr is set to zero. For boundaries where there is a convective flux out of the solution domain, the normal gradient of ePr is set to zero. For example, when the velocity is parallel to the x*-axis, the appropriate boundary conditions are epr(O, y *> = 0 ) whereas
for velocity
ePr(O,y *)
epr(x*,O)=O,
angles greater
=0,
$1
=O,
ePr(x*, 2R) = 0 ,
x*=2R
than 0” but less than 90”, the boundary
epr(X*,O)=O,
31
=o, x’=2R
%I aY
conditions
=o. y*=ZR
are
(10)
These boundary conditions have been chosen so that numerical solutions obtained using central difference schemes do not have ‘wiggles’ caused by the outflow boundary conditions as discussed by Roache [3]. The exact solution of (8) subject to these boundary conditions is given in Appendix A. The finite difference approximation of these equations, i.e., of (8) and (9) or (lo), has already been symbolically written as (5) and the solution to the finite difference equation set yields the total solution errors, ep. Since the operator error is simply the additional error caused by this finite difference approximation of 2 by L, e”pP = ep - ePr. Numerical solutions were obtained using two finite difference schemes. A linear profile in 4 between grid points for the convective and diffusive fluxes was assumed, leading to the classical central difference scheme (CDS) for one scheme. A constant stepwise profile was assumed for the convective fluxes leading to the classical upwind difference scheme (UDS) [3] for the other scheme. The exact and numerical results are presented for Peclet numbers (Pe) of 1 and 10 and for flow angles (0) of 0” and 45”. Contours of constant ePr arising from the single source have been evaluated from the exact solution and are plotted as solid lines in Figs. 3-6. Contours of constant e have been evaluated from the finite difference solutions and are plotted as broken lines. The closeness of the broken line to the solid line provides a direct indication of operator error. Since the results in each case are symmetrical about the streamline through the centre of
G.D. Stubley et al., Finite difference and influence schemes for convection and di@siun
1.58
20.07
20.0 Pe-
Ps = 1.0
1.0
8 =45*
8 = 0’
e p’=o.05
15.0 -
i
/
,
/
Y* 10.0-
5.0 -
o.oL 5.0
0.0
IO.0
15.0
0.0
20.0
5.0
I
I
I
10.0
15.0
20.0
X*
X* Fig. 3. CDS and UDS operator test for Pe = 1.0 and t? = 0”.
Fig. 4. CDS and UDS operator test for Pe = 1.0 and 0 = 45”.
the source volume only half of the contours from the numerical solutions are given. The upper halves are the contours based on CDS results and the lower halves are based on UDS results. The contours are for lines of e = 0.05. The results presented in Figs. 3 and 4 are for Pe = 1.0. For this low grid Peclet number the difference between the profile and total error (i.e., operator error) is small for both schemes. The UDS contour is slightly wider than the exact contour for 8 = 4.5”, which indicates that the influence of the residual source is being incorrectly propagated normal to the velocity vector. 20.0
i-
I
Pa = 10.0 8
= 00 IS.0
Y* 10.0
:~~,,,” “.y o.oliI 0.0
5.0
10.0
15.0
20.0
X0
Fig. 5. CDS and UDS operator test for Pe = 10.0 and i3 = 0”.
0.0
5.0
10.0
,
i
15.0
20.0
X+
Fig. 6. CDS and UDS operator test for Pe = 10.0 and 6 = 45”.
G.D. Stubley et al., Finite difference and influence schemes for convection and difision
1.59
The results are much more dramatic for the higher Peclet number, Pe = 10.0. For 0 = O”, the UDS results, shown in Fig. 5, closely match the exact contour except near the source volume. A portion of the CDS results match the exact contour downstream of the source volume but additional closed contours appear upstream of the source. This indicates that the CDS operator incorrectly propagates the influence of the residual source upstream. For 8 = 45” the upstream propagation by the CDS operator is again evident (Fig. 6); negative values of e as large as 63.5% of the value at the centre of the source volume appear. The UDS operator gives a contour that is much wider than the exact contour. This incorrect propagation normal to the velocity vector has been referred to as ‘false’ or ‘numerical diffusion’ [3, 41. This test problem permits the convection-diffusion characteristics of the CDS and UDS operators to be clearly demonstrated and understood. Other numerical operators can also be tested and compared using the same problem.
4. Performance
of influence schemes
When the finite difference equations are derived by conventional means, the choice of the profile shape fixes both the profile error and the operator error. In an attempt to assert a certain amount of separate control over each of the error components, the authors have proposed a new class of finite difference schemes called ‘influence schemes’ [2]. A different motivation led Chen and co-workers [7] to independently propose a similar class of schemes. A brief description of influence schemes and the results of subjecting two such schemes to the test problem described in the previous section will now be given. The P-point in Fig. 7 has, as its immediate neighbours, the 8 points N, NE, E, . . . , NW. An accurate simulation of the convection-diffusion process is desired within the subdomain whose boundary lies through these points (solid lines in Fig. 7) so that the +-value at each boundary point accurately ‘influences’ the value of &. In addition, the presence of sources of 4 within each of the 9 subsections (denoted by dotted lines in Fig. 7) should also accurately influence
l
BOUNDARY SUBDOMAIN /
Fig. 7. Subdomain
for determining
r-
OF
SUBDOMAIN
the influence
coefficients.
160
G.D.
Stubley et al., Finite difference and influence schemes for convection
and di#iuion
the value of &. To obtain accurate influence, the finite difference operator, L, should closely simulate the differential operator in (1). The L-operator was therefore based on an analytical solution to the equation
w [pulp ax
+
“*4+(-j [PVIP-w = rP 2a*++r dX
Pdy2
dY
(11)
r
where the coefficients are constant, evaluated at P, and d1 is constant, but perhaps different, over each of the 9 subsections in Fig. 7. To solve (11) continuous boundary conditions around the subdomain in Fig. 7 are required. This requires a profile assumption to specify a distribution of 4 between the grid points around the subdomain. Linear and quadratic profiles give rise respectively to the linear influence scheme (LIS) and the quadratic influence scheme (01s). The resulting finite difference equation for the P-subdomain has the form
where the & are the values of 4 at the 8 neighbour points, and the the 9 subsections. The analytical solution of (11) subject to the ditions, results in expressions for the ‘influence coefficients’, C1, and coefficients’, DI. The derivation of these expressions, by the Green’s underlines the influence point of view [6]. The expressions for Appendix B. When the two influence schemes, LIS and QIS, are applied to the the last section, with Pe = 1.0 and 13= 0” and 45”, the numerical and
0, are the source terms in prescribed boundary conthe ‘source term influence function method, further C, and Dr are given in test problem, described in analytical solutions (i.e., e
20.0 -
20.0-
Pe = 10.0
PO = 10.0
15.0 -
a. I.S.--eP’=0.05
Y* 10.0-
0.0 0.0
I 5.0
I 10.0
I 15.0
I 20.0
X*
Fig. 8. QIS and LIS operator test for Pe = 10.0 and 6 = 0”.
0.0 / 0.0
I 5.0
I 10.0
1 15.0
I 20.0
X*
Fig. 9. QIS and LIS operator test for Pe = 10.0 and 0 = 45”.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
161
and epr) were virtually identical. For Pe = 10.0, Figs. 8 and 9 show that the QIS leads to for 8 = 0” and 45”, respectively. These figures also show that the LIS accurate predictions operator gives a contour slightly wider than the exact contour for 8 = 45”. It should also be noted that small negative values (3.4% of the value at the centre of the source volume) appear in the QIS solution for Pe = 10.0 and 8 = 45”. These results clearly indicate that the operator errors of the influence schemes have been kept to a low level. While the operator error is principally determined by the operator in (11) the results indicate that the profile approximations still play a slight role in determining the operator error.
5. Improvements
to the influence schemes
Two approximations are required to derive the finite difference equations for the influence schemes, a profile approximation around the subdomain (Fig. 7) and the form of the operator used to obtain the influence coefficients, (11). In the previous section and earlier work [2], a comparison of results based on linear and quadratic profiles has shown that improving the profile approximation reduces solution error. In this section it will be shown that the finite difference operator derived from (11) is still inadequate for some problems, and that a better operator can easily be derived. The finite difference operator for the influence schemes, and Chen’s scheme, is based on an analytical solution of (11) as described. The constant coefficients, assumed in (11) can lead to large errors when the velocity and/or diffusion coefficients vary sharply over the subdomain (Fig. 7). This error can be substantially reduced by modifying the source terms, as will now be shown. Eq. (1) which defines the exact operator, can be written without approximation as
[pfi]p
g
+
[pv’],
a*@ +r”*@ E = r---y 7+d ax dY aY
where
o= d-Cori-[pa],)~-@6-[p6]~)~.
(W
The mass conservation equation has been used to modify the convection terms, the diffusion terms have been expanded, and new apparent mass fluxes have been defined as pti = PU - arlax,
pc = PV - arjay
.
(14)
Eq. (13a) is of the same form as (11) if F is approximated by r, and if 0 is taken as constant over each of the 9 subsections in Fig. 7. Thus expressions for C, and D1 in (12) based on (13a) are identical to those based on (11) except that [piilP and [p& replace [pulp and [pv)lP, respectively. In evaluating d for each subsection, the values of pfi, pv’, a@/dx and a@/ay are evaluated at the grid point associated with that subsection; &?/ax and a@/ay are based on the best available (usually from the last iteration) values of 4 and the derivatives are ap-
162
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
proximated by passing quadratic profiles through the surrounding grid points. The linear and quadratic influence schemes that have been extended in this manner have been named LISE and QISE respectively. The importance of these extensions can be seen by applying the influence schemes to obtain numerical solutions to the following problem, proposed by Baliga and Patankar [8]. The distribution of a conserved scalar @ in a flow between a rotating and fixed cylinder, Fig. 10, is sought on a Cartesian grid within the test section shown in the north-east quadrant. In Cartesian coordinates the problem appears to be two-dimensional with velocity and diffusivity fields which may vary spatially. Analytical expressions for the distributions of the velocity and the scalar @ can be readily obtained by formulating the problem in cylindrical coordinates. The conditions under which solutions were obtained in the present work are: (i) the outer cylinder is fixed and the inner cylinder rotates at an angular velocity wi; (ii) all variables are nondimensionalized by their values at the inner cylinder radius ri, i.e.,
(iii) the diffusivity of the scalar @ is given by r* = 1.0/r*2; (iv) @ is 1.0 at the inner cylinder and 0.0 at the outer cylinder. In cylinder coordinates the expression for the velocity is V,*(r*) =
(r*2 - ro*‘) r*(l.O - rZ2)
and that for the scalar @ is
&.‘I,; Fig. 10. Layout
for rotating
the cylinder
test problem.
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
163
(17) For this problem
(1) simplifies
to
where
U
*=
“*yT f3 f%
3
V
*-&
pe=
”
(19)
pirToi’
Numerical solutions to (18) were obtained using the quadratic influence scheme, QIS, the extended influence scheme, QISE, and the exponential difference scheme, EDS [2], which reduces to CDS at low Peclet numbers and to UDS at high Peclet numbers. The exact solution, (17) was used to specify the values on the boundary of the test section. The outer cylinder radius, r i, was set at 3.0 and each scheme was tested for a wide range of Peclet number. The set of finite difference equations, for each test, was solved using the line-by-line method [l], and the solution was considered to have converged when the maximum change of 4 at every grid point was less than 0.001% between iterations. The performance of each scheme was assessed by studying the percentage error in @, that was predicted at the centre of the test section, as the Peclet number was varied. The results, Fig. 11, demonstrate that the solution based on the influence scheme, QIS (without modification), had prohibitively large errors for this problem at small Pe (when diffusive transport dominates). The extended influence scheme, QISE, outperforms, in terms of
VI0
=
0.0
r = l/r’*
0.01
0.1
IO
1.0
100
1000
Pe+ Fig. 11. Midpoint
errors
for the rotating
cylinder
test problem.
164
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
accuracy, both EDS and QIS for all Peclet numbers tested. Similar results were obtained for this problem but under different conditions i.e., where both cylinders rotated at the same angular velocity, or the diffusivity was constant, etc. [6]. These results show the large solutions errors that can arise by using QIS in cases where the velocity and diffusivity vary spatially. By using the QISE, the solution errors can be greatly reduced with little additional numerical effort.
6. Summary This paper makes two separate contributions and, because the basic ideas will be new to most potential readers, a review of some previous work was also presented. First the equations for the error in a finite difference solution were written and the motivation for subdividing the error into ‘profile’ and ‘operator’ components was reviewed. The truncation error of the finite difference scheme provides a qualitative measure of the profile error but no good test for operator error was available in the literature. The first main contribution of this paper was to propose such a test and to subject a central difference scheme, an upwind difference scheme, and two influence schemes to this test; the ideas behind the formulation of influence schemes were reviewed. While influence schemes have desirable properties, it is shown that large solution errors can still exist if the diffusion coefficients and velocities change rapidly over the solution domain. The second contribution of this paper is to show a simple procedure for greatly reducing the solution error for such problems.
Acknowledgment The financial support of the Natural Sciences and Engineering Research Council of Canada, through operating grants and through a scholarship to G.D. Stubley, is gratefully acknowledged.
Appendix A In the text a test for the operator error was discussed. The profile error, solving (8) subject to the boundary conditions given in (9) or (lo), depending the velocity vector. The solution to (8) given in this appendix is obtained function method [6]. The exact solution to (8) subject to the boundary conditions given in (9)
epyx*,y*) = 2
sir@y *) exp[Pe,x
“=l y R + $$ where
*/2 + Pe,y */2] I&
sin(2@?))(cosh(2yR)
+2
sinh(2yR))
epr, is obtained by upon the angle of using the Green’s is
(A4
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
p = nn/2R,
y2 = Pe2/4 + p” ,
Pe,=Pecosf3,
Pe, = Pe sin
165
8,
Rttl2 II
=
I
R-112
evPe,d21 sin(h) dv , sinh(yx*)[cosh(y(2R
- 5)) + %
sinh(y(2R
- &))I ; 6 > **
- x*)) + %
sinh(y(2R
- x*))]
expPe,5/21
d5.
sinh(y5)[cosh(y(2R
;e < x*
I
The exact solution to (8) subject to the boundary conditions given in (10) is identical given in (A.l) except for the eigenvalues which are now given as p = (2n - 1)7r/4R.
Appendix
to that
(A.21
B
The influence coefficients, been obtained [6] by solving expressions are given in this The influence coefficients r = AxlAy
C,..,, etc., and the source term influence coefficients, D,, etc., have (11) analytically using the Green’s function method. The resulting appendix. are functions of three nondimensional parameters:
(grid aspect
ratio),
Pe, = [PU],AX/T~
(local x Peclet
number)
Pe, = [p~]~Ay/r~
(local y Peclet
number).
The resulting
expressions
for each influence -
ci =exdipex+ iPeylmsl 2
coefficient
(-l)“Af,, cosh(a r) m
(B-1)
,
have the form
03.2)
where u’, = a2 + pz, &, = i(2m - 1)~ and a2 = Pez/4r2 + Pet/4. The expressions for the series coefficients, A”, and A”,“, are given in Table B. 1 for LIS and QIS. It should be noted that influence of each corner point is broken into two components so that
G’ = Gvw +
csws
(B.3)
where the expressions for C,, and C,,, have the form given in (B.2). The expressions for the is other influence coefficients are all related to C, and C,,,. A summary of these relationships given in Table B.2. The expressions for the source term influence coefficients have the form
G.D. Stubley et al., Finite difference and influence schemes for convection and diffusion
166 Table B.1 Expressions
for the series coefficients
Influence coefficient
Boundary condition function
CW
f = {(2’;:,,
o
Am for the influence
Scheme 1
schemes
Series coefficient
LIS where
y = :Pe: + /3$
QIS
f=Y*(2-y*),
CW
Ah
O
C SWW
1
C SW
h =(l-y*)(2-y*)/2,
LIS QIS
k-3+_
Y
peP+&($peLP2) Y Y3
m
o
03.4) where the series coefficient pairs, I&, for each influence coefficient are
The series coefficients, Ii and lj, are given in Table B.3. To aid in the convergence of the series (B.4) for certain coefficients the following identities have been derived.
&
D,
r 0, AY2
=epey"S2 - e Pei4(iPeyS0 + S,) + exp[iPe,
+ tPe,] .z, 2J-~~~$$ m
), m
= eP’~‘4($Pe,(1 - e-Pey’2)So + (1 + e-pey’z)S1)
I, &
D, = e-Pey’2S2+ epPe~4(~Pe,,So- S,) + exp[$Pe, + $Pe,] 2 (- ‘)“lrxr, m=, 2~~ cos(ra,)
03.6)
’
167
G.D. Stubley et al., Finite difference and influence schemes for conuection and diffusion
where 2 SO= Pe$ sinh(Pe,)
coth(Pe,) + &-)(cosh(sPe,)
- cosh(iPe,))
Y
- %sinh(jPe,) + j sinh(iPe,) ,
s, = Pe
siih(Pe Y
) coth(Pe,)(sinh(!Pe,)
+ sinh($r
Y
- i sinh($Pe,) - a cosh(iPe,) , 1 sinh(PeJ2) ” = 2Pe, cosh’(PeJ2)
(3.7)
’
Table B.2 Summary of influence coefficients in terms of functions for C, and C,, C, = Csw(Pex, Pe,, r) + CswPe,, Pex, l/r) -Pe,, l/r) + e-“xC,(PeX. Pe,, r) C,, = Cm(Pex, -Pe,, r) + ewpeyG.&k,, Pe,, I/r) C,, = e-P”YCsww(Pey,-Pe,, l/r) + e-PeW~w(Per, -Pe,,
C, = C,(Pe,, Pe,, r)
c,= Csw(PeY,
Table B.3 Expressions
C, = e-P’xCw(Pe,, Pe,, r) C, = C,(Pe,, -Pe,, l/r)
r)
C, = ewPeyCw(Pey,-Pe,,
for source term influence coefficients
= 2u,,, cosh(ra,)
I PX
e-~ex~2
+
I n
Y
e -3Pe&T
I, = - r/ - c k
sinh $rcr,,,- ra, co& k,,,]
[$e, [$e,
e-3PQl‘i = ry
[-$e,
sinh dru,,, + ru,,, cash sinh
irm,,,
+
?ra,]
ra, cash ircr,,,] - a, Y
&+
4v
=
COS(&~~)[&
F ewPey4 CoS(b> Y
+ iPe,(-
[1p (_I)“(1 2 ey
_
I)“1
e-Pey/Z)
+ pm(l
+ e-*W)]
l/r)
168
G.D. Stubley et al., Finite difference and influence schemes for convection
and diffusion
References [l] S.V. Patankar, Numerical Heat Transfer and Fluid Flow (McGraw-Hill, New York, 1980). [2] G.D. Stubley, G.D. Raithby and A.B. Strong, Proposal for a new discrete method based on an assessment of discretization errors, Numer. Heat Transfer 3 (1980) 411-428. [3] P.J. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, NM, 1972). [4] G.D. Raithby, A critical evaluation of upstream differencing applied to problems involving fluid flow, Comput. Meths. Appl. Mech. Engrg. 9 (1976) 7.5-103. [5] K. Torrance et al., Cavity flows driven by buoyancy and shear, Fluid Mech. 51 (1972) 221-231. [6] G.D. Stubley, A new discrete method for convection-dominated flows based on a clear understanding of solution errors, Ph.D. Thesis, University of Waterloo, Waterloo, Ont., Canada, 1981. [7] C.J. Chen and P. Li, Finite differential method in heat conduction-application of analytical solution technique, ASME Paper 79-WA/HT-50, 1979. [8] B.R. Baliga and S.V. Patankar, A new finite-element formulation for convection-diffusion problems, Numer. Heat Transfer 3 (1980) 395410.