����������� ������ 47 (1995), 103–111

UDC 519.633

���������� ������ ��� research paper

FLOW-ORIENTED DIFFERENCE SCHEME FOR MULTIDIMENSIONAL CONVECTION-DIFFUSION EQUATION

Desanka Radunović

Abstract. In this paper one flow-oriented difference scheme for multidimensional convection-diffusion equation is constructed and analysed. The order of the accuracy is O(Δx^{2}), except for convection dominant case when it decreased by one. The stability depends on the diffusion coefficient D, and for the square grid the stability condition is D Δt/Δx^{2} ≤ 0.25. Some examples are presented to illustrate that the scheme is especially applicable for dominantly convection problems and problems with not enough smooth solutions.

1. Introduction

The convection-diffusion equation is one of the fundamental equations of microelectronics. In the last several years there has been considerable effort aimed at semiconductor device modelling ([1],[2]). Standard finite-difference or finite element methods give rather poor results in modelling such problems, especially when the convection term is dominant, because spurious oscillations appear. The equation is officially parabolic, but practically hiperbolic as the diffusion term is negligible in comparison with the convection one. Thus, interior and boundary layers, i.e. thin regions of fast variation of the solution, appear. This causes difficulties which can hardly be surpassed by standard techniques. In [3] the necessary criteria which a good numerical discretization scheme for the convection-diffusion equation must satisfy are identified. Various ideas were implemented to develop stable and accurate methods for solving these problems ([3]-[6]).

In this paper one flow-oriented difference scheme for multidimensional convection-diffusion equation is constructed and analysed (similar approach can be found in [7]). The basic idea is given in Section 2 and it is applied and analysed for multidimensional problem in Section 3. The order of the accuracy is O(Δx^{2}), except for convection dominant case when it decreased by one. The stability criterion is derived in Section 4. In the last section some examples are presented to illustrate that the scheme is especially applicable for dominantly convection problems and problems with not enough smooth solutions.

1991 AMS Subject Classification: 65C20, 76R99

2. Onedimensional problem

Let us begin with the onedimensional problem

\frac{d}{dx}(D \frac{du}{dx} − Cu) = R,  x ∈ (0, 1),  (1)

D, C and R generally depends on x and u, with arbitrary boundary conditions. The problem (1) is discretised on the uniform grid

ω_{x} = \{x_{i} | x_{i} = iΔx, Δx = \frac{1}{n}, i = 0,…, n\}.

By use of notation

V_{i+ \frac{1}{2}} = (D \frac{dv}{dx} − C v)|_{x=x_{i+ \frac{1}{2}}} ,  (2)

where v(x) is the approximate solution of (1), and

δ_{x}f(x) = \frac{1}{Δx} (f(x + \frac{Δx}{2}) + f(x − \frac{Δx}{2})) ,

the conservative difference approximation of the equation (1) in the interior node x_{i} is

δ_{x}V_{i} = R(x_{i}, v(x_{i})).  (3)

If the coefficients D and C of the equation (1) are approximated part by part by constants,

D(x) ≡ D_{i+ \frac{1}{2}} = D(x_{i+ \frac{1}{2}} , u(x_{i+ \frac{1}{2}} )), C(x) ≡ C_{i+ \frac{1}{2}} = C(x_{i+ \frac{1}{2}} , u(x_{i+ \frac{1}{2}} )),

x ∈ (x_{i}, x_{i+1}), the relation (2) for x ∈ (x_{i}, x_{i+1}) can be treated as the first-order differential equation with constant coefficients, and its solution is

v(x) = exp (\frac{C_{i+ \frac{1}{2}}}{D_{i+ \frac{1}{2}}} (x − x_{i+ \frac{1}{2}} )) ×

(v(x_{i+ \frac{1}{2}} ) + \frac{V_{i+ \frac{1}{2}}}{C_{i+ \frac{1}{2}}} − \frac{V_{i+ \frac{1}{2}}}{C_{i+ \frac{1}{2}}} exp ( − \frac{C_{i+ \frac{1}{2}}}{D_{i+ \frac{1}{2}}} (x − x_{i+ \frac{1}{2}} ))) . (4)

By eliminating v(x_{i+ \frac{1}{2}} ) from (4) by use of v(x_{i}) and v(x_{i+1}) we can express the term which we need for the difference approximation (3)

V_{i+ \frac{1}{2}} = −\frac{2}{Δx}D_{i+ \frac{1}{2}} α_{i+ \frac{1}{2}} \frac{v_{i+1} exp(−α_{i+ \frac{1}{2}} ) − v_{i} exp(α_{i+ \frac{1}{2}} )}{exp(−α_{i+ \frac{1}{2}} ) − exp(α_{i+ \frac{1}{2}} )}.  (5)

Here, v_{i} = v(x_{i}) and α_{i+ \frac{1}{2}} = Δx C_{i+ \frac{1}{2}} /(2 D_{i+ \frac{1}{2}} ) is the mesh Peclet number, which controls the rate of mixing on the grid element. Thus from (5) and (3) the conservative difference scheme for the equation (1) becomes

\frac{D_{i+ \frac{1}{2}}}{Δx^{2}} (α_{i+ \frac{1}{2}} cth(α_{i+ \frac{1}{2}} )(v_{i+1} − v_{i}) − α_{i− \frac{1}{2}} cth(α_{i− \frac{1}{2}} )(v_{i} − v_{i−1}))

− \frac{1}{2Δx} (C_{i+ \frac{1}{2}} (v_{i+1} + v_{i}) − C_{i− \frac{1}{2}} (v_{i} + v_{i−1})) = R_{i}. (6)

If the coefficients D and C are constant and R ≡ 0, the solution of the difference equation (6) v_{i} = A+B exp(iΔx C/D) (A and B are appropriate constants) agreeds in the grid nodes with the solution of the equation (1) u(x) = A + B exp(x C/D).

For the constant coefficients central difference scheme

δ_{x}(Dδ_{x}v − C μ_{x}v)|_{x=x_{i}} = 0,  (7)

where μ_{x}v(x) = \frac{1}{2} (v(x + \frac{Δx}{2}) + v(x − \frac{Δx}{2})), the solution is

v_{i} = A + B(1 + \frac{Δx}{2} \frac{C}{D} )^{i}(1 − \frac{Δx}{2} \frac{C}{D} )^{−i}.

Thus, it is obvious that for singularly perturbed problems (small D) the solution is unstable if the mesh size is not small enough.

The expression (5) can be rewritten as

V_{i+ \frac{1}{2}} = (D δ_{x}v + D_{α} δ_{x}v − C μ_{x}v)|_{x= \frac{1}{2} (x_{i}+x_{i+1})} ,  (8)

where D_{α} = D(α_{i+ \frac{1}{2}} cth(α_{i+ \frac{1}{2}} )−1). By comparison of the schemes (7) and (3),(8), it is obvious that the diffusion correction D_{α} δ_{x}v made the second scheme stable also for the large Peclet number α, i.e. when the convection is dominant.

3. Multidimensional problem

Now, consider the multidimensional convection-diffusion equation. Direct implementation of the scheme (6) to every coordinate direction do not give satisfactory results, as a nonphysical diffusion in the direction normal to the convection, so called ”crosswind“ diffusion, produces spurious oscillations. One stable and nondissipative scheme is obtained by use in the convection direction the diffusion correction described in the previous section.

Without lose of generality, we shall treat the twodimensional problem

Lu ≡ \frac{∂}{∂x} (D \frac{∂u}{∂x} − C_{x}u) + \frac{∂}{∂y} (D \frac{∂u}{∂y} − C_{y}u) = R  (9)

subject to boundary conditions, on the unit square. We construct the difference scheme for the equation (9) on the rectangular grid

ω_{xy} = \{(x_{i}, y_{j}) | x_{i} = iΔx, y_{j} = jΔy, Δx = \frac{1}{n}, Δy = \frac{1}{m}, i = 0,…, n, j = 0,…, m\},

To obtain the diffusion correction to the convection direction we rotate the coordinate system (x, y) for the angle

θ = arctg (C_{y}(x_{i}, y_{j}, u(x_{i}, y_{j}))/C_{x}(x_{i}, y_{j}, u(x_{i}, y_{j}))) ,
[ \begin{matrix}ξ \\ η \end{matrix} ] = [ \begin{matrix}cos θ & sin θ \\ − sin θ & cos θ \end{matrix} ] [ \begin{matrix}x \\ y \end{matrix} ] ,  (10)

for every grid node (i, j) locally. In new coordinates, assuming that the coefficients of the equation are constant, the equation (9) is

D (\frac{∂^{2}u}{∂ξ^{2}} + \frac{∂^{2}u}{∂η^{2}}) − \frac{∂}{∂ξ} (Cu) = R, (11)

where

C = \sqrt{C_{x}^{2} + C_{y}^{2}} sign(C_{x}C_{y}).

To apply the approximation (6) we shall rewrite (11) as

D (\frac{∂^{2}u}{∂ξ^{2}} + \frac{∂^{2}u}{∂η^{2}}) − D \frac{∂^{2}u}{∂ξ^{2}} + \frac{∂}{∂ξ} (D \frac{∂u}{∂ξ} − Cu) = R.  (12)

The Laplace operator is invariant towards the rotation (10), so is

D (\frac{∂^{2}u}{∂ξ^{2}} + \frac{∂^{2}u}{∂η^{2}}) ≡ D (\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}) ≈ D (δ_{x}^{2}v + δ_{y}^{2}v) .  (13)

Only derivatives to the ξ direction appear in the remaining two terms of the left side of the equation (12), thus the second term can be approximated by the onedimensional central difference

D \frac{∂^{2}u}{∂ξ^{2}} ≈ Dδ_{ξ}^{2}v ≡ \frac{D}{ΔL^{2}} (v_{D} − 2v_{i,j} + v_{U}),  (14)

and the last one by use of the appropriate form of the scheme (6). Approximate values of the solution in the points D and U (fig.1.) are determined by the bilinear interpolation from four the nearest nodes of the grid ω_{xy}:

\begin{matrix}v_{D} & = v_{i+p,j} + \frac{Δx}{Δy}(v_{i+p,j+q} − v_{i+p,j})| tg θ| \\ v_{U} & = v_{i−p,j} + \frac{Δx}{Δy}(v_{i−p,j−q} − v_{i−p,j})| tg θ| \\ ΔL & = Δx\sqrt{1 + tg^{2} θ} \end{matrix} | tg θ| ≤ \frac{Δy}{Δx},

or

\begin{matrix}v_{D} & = v_{i,j+q} + \frac{Δy}{Δx}(v_{i+p,j+q} − v_{i,j+q})| ctg θ| \\ v_{U} & = v_{i,j−q} + \frac{Δy}{Δx}(v_{i−p,j−q} − v_{i,j−q})| ctg θ| \\ ΔL & = Δy\sqrt{1 + ctg^{2} θ} \end{matrix} | tg θ| > \frac{Δy}{Δx},

where p = sign(C_{x}) and q = sign(C_{y}).

We replace (13), (14) and (6) (where we put v_{i+1} ≡ v_{D}, v_{i−1} ≡ v_{U}) into (12) and obtain the difference approximation of the eq.(11) in the interior node of the grid ω_{xy}

Λv ≡ D (δ_{x}^{2}v + δ_{y}^{2}v) + D (α cth(α) − 1) δ_{ξ}^{2}v − Cμ_{ξ}v,  (15)

where α = C ΔL/(2D).

For constant coefficient case and smooth boundary conditions the consistency of the approximation (15) follow from the expansion

Lu − Λu = D(\frac{Δx}{ΔL})^{2}(α cth(α) − 1) (\frac{∂^{2}u}{∂x^{2}} + 2\frac{C_{y}}{C_{x}} \frac{∂^{2}u}{∂x∂y} + \frac{C_{y}}{C_{x}} \frac{Δy}{Δx} \frac{∂^{2}u}{∂y^{2}}) + O(Δx^{2})
= const ∗ \frac{Δx^{2}}{D} + O(Δx^{2}) = \{ \begin{matrix}O(Δx), & for D = O(Δx), \\ O(Δx^{2}), & for D = O(1). \end{matrix}  (16)

As α cth α − 1 = O(α^{2}) = O((ΔL/D)^{2}) it is clear that for singularly perturbed problems (for small D) the order of the accuracy is decreased by one.

4. Stability of the initial-value problem

We approximate the nonstationary convection-diffusion equation

\frac{∂u}{∂t} − Lu = R,  (17)

where the operator L is defined in (9), by the second-order accurate approximation (15) accomplished via the Euler forward scheme for differentiation in time

δ_{t}v = Λv + R.  (18)

We put the Fourier expansion of the discrete solution

v_{i,j}^{n} = ∑_{m_{1},m_{2}} ζ^{n} exp(ı(m_{1} i Δx + m_{2} j Δy)),  (19)

to the difference equation (18) and obtain that the amplifying factor is

ζ =1 − 4σ(\frac{C_{x}^{2}}{C_{x}^{2} + C_{y}^{2}} (α cth(α) − 1)(K sin^{2} \frac{m_{1}Δx + m_{2}Δy}{2}

+ (1 − K) sin^{2} \frac{m_{1}Δx}{2}) + sin^{2} \frac{m_{1}Δx}{2} + (\frac{Δx}{Δy})^{2} sin^{2} \frac{m_{2}Δy}{2} )

− ıC_{x}\sqrt{\frac{σΔt}{D}}(K sin(m_{1}Δx + m_{2}Δy) + (1 − K) sin(m_{1}Δx)),

where σ = D Δt/Δx^{2} and K = Δx tg(θ)/Δy. Therefore, we can write

|ζ|^{2} = f_{0}(m_{1}Δx, m_{2}Δy) + Δt f_{1}(m_{1}Δx, m_{2}Δy) + Δt^{2} f_{2}(m_{1}Δx, m_{2}Δy),

where

f_{0}(m_{1}Δx,m_{2}Δy) = (1 − 4σ( sin^{2} \frac{m_{1}Δx}{2} + (\frac{Δx}{Δy})^{2} sin^{2} \frac{m_{2}Δy}{2} ))^{2},
f_{1}(m_{1}Δx,m_{2}Δy) = 8\frac{C_{x}^{2}}{C_{x}^{2} + C_{y}^{2}} \frac{σ(α cth(α) − 1)}{Δt} (K sin^{2} \frac{m_{1}Δx + m_{2}Δy}{2}
+(1 − K) sin^{2} \frac{m_{1}Δx}{2}) (4σ( sin^{2} \frac{m_{1}Δx}{2} + (\frac{Δx}{Δy})^{2} sin^{2} \frac{m_{2}Δy}{2} ) − 1)

+ C_{x}^{2} \frac{σ}{D} (K sin(m_{1}Δx + m_{2}Δy) + (1 − K) sin(m_{1}Δx))^{2},

f_{2}(m_{1}Δx,m_{2}Δy) = (4\frac{C_{x}^{2}}{C_{x}^{2} + C_{y}^{2}} \frac{σ(α cth(α) − 1)}{Δt}(K sin^{2} \frac{m_{1}Δx + m_{2}Δy}{2}

+(1 − K) sin^{2} \frac{m_{1}Δx}{2}))^{2} .

f_{1}(m_{1}Δx, m_{2}Δy) and f_{2}(m_{1}Δx, m_{2}Δy) are bounded functions, thus the approximation (18) will be stable for |f_{0}(m_{1}Δx, m_{2}Δy)| < 1, i.e. for

DΔt ≤ \frac{1}{2} \frac{(ΔxΔy)^{2}}{Δx^{2} + Δy^{2}} ,  (20)

or, for the square grid (Δx = Δy), for D \frac{Δt}{Δx^{2}} ≤ \frac{1}{4}.

The convergence of the scheme (18) follows from the Lax’s equivalence theorem ([8]).

5. Numerical results

Following examples illustrate stability and accuracy of the flow-oriented scheme (FOS) for linear and also for nonlinear problems. For the comparison results calculated by use of the central finite difference scheme (CD),

δ_{t}v_{i,j}^{n} = D(δ_{x}^{2}v_{i,j}^{n} + δ_{y}^{2}v_{i,j}^{n}) + C_{x} δ_{x}v_{i,j}^{n} + C_{y} δ_{y}v_{i,j}^{n},

and the upwind scheme (UP)

δ_{t}v_{i,j}^{n} = D(δ_{x}^{2}v_{i,j}^{n} + δ_{y}^{2}v_{i,j}^{n}) + \frac{C_{x}}{2} ((1 − p)δ_{x}v_{i+ \frac{1}{2} ,j}^{n} + (1 + p)δ_{x}v_{i− \frac{1}{2} ,j}^{n})

+ \frac{C_{y}}{2} ((1 − q)δ_{y}v_{i,j+ \frac{1}{2}}^{n} + (1 + q)δ_{y}v_{i,j− \frac{1}{2}}^{n} )

where p = sign(C_{x}) and q = sign(C_{y}), are presented. In all examples the square grid is used.

Example 1. The problem is defined on the unit square by the equation

\frac{∂u}{∂t} = D (\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}) − cos(θ) \frac{∂u}{∂x} − sin(θ) \frac{∂u}{∂y}.

The initial conditions are u = 1 for 0 ≤ x ≤ 1, 0 ≤ y ≤ (1 − x)/4, otherwise u = 0. The boundary conditions are defined by the initial one and do not change in time.

Calculations are done for grid parameters Δx = Δy = 0.05 and time step Δt = 0.002. The direction of the convection is determined by the choise of the angle θ. FOS scheme diverges for D = 10, that agreed with (20) (D Δt/Δx^{2} = 8), and converges for small D. The results for θ = 45^{∘} and D = 10^{−3} obtained by use of FOS (fig.2), CD (fig.3.) and UP (fig.4.) schemes are given. The dissipation of UP scheme is obvious, while CD scheme is very unstable. The nonlinear filtering technique proposed in [9] is applied to stabilize CD scheme (fig. 5.), but the results obtained by FOS scheme are still much better.

Fig. 2. FOS scheme Fig. 3. CD scheme

Fig. 4. UP scheme Fig. 5. CD + filter scheme

Example 2. The linear problem, defined on the unit square by the equation

\frac{∂u}{∂t} = D (\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}) − w(x, t) \frac{∂u}{∂x} − w(y, t) \frac{∂u}{∂y},

where w(x, t) = (0.1A + 0.5B + C)/(A + B + C) and A = exp(−0.05(x − 0.5 + 4.95t)/D), B = exp(−0.25(x−0.5+0.75t)/D) and C = exp(−0.5(x−0.375)/D), is solved for t ∈ (0, 1) and D = 10^{−4}. The initial and Dirichlet boundary conditions are determined by the true solution u(x, y, t) = w(x, t) w(y, t).

Fig. 6. presents FOS solution of the problem for t = 0.6 calculated for parameters Δx = Δy = 0.05 and Δt = 0.002. The comparison of the numerical solution v(x, y, t) for y = 0.6 and t = 0.6 obtained for various steps Δx = Δy and constant ratio Δt/Δx^{2} = 0.8 is given on fig.7. The convergence is obvious.

Fig. 6. FOS solution for t = 0.6 Fig. 7. FOS solution for y = 0.6 and t = 0.6

Example 3. The nonlinear problem, defined on the unit square by the equation

\frac{∂u}{∂t} = D (\frac{∂^{2}u}{∂x^{2}} + \frac{∂^{2}u}{∂y^{2}}) − 3u \frac{∂u}{∂x} − 3(1.5 − u) \frac{∂u}{∂y},

is solved for t ∈ (0, 1) and D = 3 ∗ 10^{−4}. The initial and Dirichlet boundary conditions are determined by the true solution

u(x, y, t) = \frac{1}{4} (3 − \frac{1}{1 + exp(0.375(−x + y − 0.75t)/D)}) .

Fig. 8. FOS solution for t = 0.6 Fig. 9. FOS solution for y = 0.6 and t = 0.6

Fig. 8. presents FOS solution of the problem for t = 0.6 calculated for parameters Δx = Δy = 0.05 and Δt = 0.002. The comparison of the numerical solution v(x, y, t) for y = 0.6 and t = 0.6 obtained for various steps Δx = Δy and constant ratio Δt/Δx^{2} = 0.8 is given on fig. 9. The delay of the numerical solution irrespective of the mesh size, not perceptible for the linear problem (fig. 6.), is obvious for this nonlinear problem. This might be caused by the low accurate explicit difference scheme used for time differentiation.

REFERENCES

[1] Mock, M.S., Analysis of mathematical models of semiconductor devices, Boole Press, Dublin, 1983.

[2] Polak, S.J., Den Heijer, C., Schilders, W.H.A., Markowich, P., Semiconductor device modelling from the numerical point of view, Int. J. Numer. Meth. Eng. 24 (1987), 763–838.

[3] Patel, M.K., Cross, M., Markatos, N.C., An assessment of flow oriented schemes for reducing false diffusion, Int. J. Num. Meth. Engrg. 26 (1988), 2279–2304.

[4] Brooks, A.H.,Hughes, J.R., Streamline upwind Petrov Galerkin formulations for convection dominated flow with particular emphasis on the incompressible Navier-Stocks equation, Comput. Meth. Appl. Mech. Engrg. 32 (1982), 199–212.

[5] Leonard, B.P., Simple high-accuracy resolution program for convective modelling of discontinuities, Int. J. Numer. Meth. Fluids 8 (1988), 1291–1318.

[6] Rai, M.M., Moin, P., Direct simulations of turbulent flow using finite-difference schemes, J. Comp. Phys. 96 (1991), 15–53.

[7] Yie, H., Zhimeng, T., Nonoscillatory streamline upwind formulations for drift-diffusion equation, IEEE Trans. Computer-Aided Design 12 (1993), 1535–1541.

[8] Richtmyer, R.D., Morton, K.W., Difference methods for initial-value problems, Interscience publishers, New York, 1967.

[9] Engquist, B., Lötstedt, P., Sjögreen, B., Nonlinear filters for efficient shock computation, Math. of Comp., 52 (1989), 509–537.

(received 27.11.1995.)

Faculty of Mathematics, Studentski trg 16, P.O.Box 550, 11000 Beograd, Yugoslavia