Unconditionally Monotone Schemes for Unsteady ...

Computational Methods in Applied Mathematics Vol. 13 (2013), No. 2, pp. 185–205 c 2013 Institute of Mathematics, NAS of Belarus

Doi: 10.1515/cmam-2013-0002

Unconditionally Monotone Schemes for Unsteady Convection-Diffusion Problems Nadyezhda M. Afanas’eva · Alexander G. Churbanov · Petr N. Vabishchevich

Abstract — This paper deals with constructing monotone schemes of the second-order accuracy in space for transient convection-diffusion problems. They are based on a reformulation of the convective and diffusive transport terms using the convective terms in the divergent and nondivergent forms. The stability of the difference schemes is established in the uniform and L1 norm. For 2D problems, unconditionally monotone schemes of splitting with respect to spatial variables are developed. Unconditionally stable schemes for problems of convection-diffusion-reaction are proposed, too. 2010 Mathematical subject classification: 65M06, 65M08, 76M12, 76M20, 76S05. Keywords: Convection-Diffusion Problems, Finite Difference Schemes, Logarithmic Norm, Monotone Schemes, Splitting Schemes.

Introduction Convection-diffusion problems are typical for mathematical models of fluid mechanics. Heat transfer as well as impurities spreading occur not only due to diffusion, but result also from medium motion. Principal features of physical and chemical phenomena observing in fluids and gases [2,11] are generated by media motion resulting from various forces. Computational algorithms for the numerical solution of such problems are of great importance; they are discussed in many publications (see, e.g., [8, 15]). In considering the second-order parabolic equations, the convective terms may be written in divergent, nondivergent, and skew-symmetric forms [23]. Transient problems of convection-diffusion are governed by evolutionary operator equations in the appropriate spaces. Their study is based on examining the corresponding properties of the differential operators of convective and diffusive transport. In constructing their discrete analogs, we focus on the approximations that inherit the basic properties of these operators. Considering boundary value problems for elliptic and parabolic equations of second order, emphasis is on the maximum principle [16]. Using this principle, there were constructed a priori estimates for convection-diffusion equations in L∞ (Ω), where the convective terms Nadyezhda M. Afanas’eva North-Eastern Federal University, 58, Belinskogo, 677000 Yakutsk, Russia E-mail: [email protected]. Alexander G. Churbanov Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tul’skaya, 115191 Moscow, Russia E-mail: [email protected]. Petr N. Vabishchevich Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tul’skaya, 115191 Moscow, Russia E-mail: [email protected]. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

186

Nadyezhda M. Afanas’eva et al.

were written in nondivergent form [10]. Similarly, in [23], the maximum principle has been formulated for problems with convective terms in the divergent form. In this case, the relevant a priori estimates hold in the space L1 (Ω). Difference schemes that satisfy the maximum principle are referred to as monotone schemes. It is well known that the simplest way to construct unconditionally monotone schemes is based on using upwind differences for the approximation of convective terms [19]. The main drawback of such schemes is associated with the first-order approximation in space. Various approaches are employed in CFD for developing monotone second-order schemes [26, 29]. Basically, the idea of switching the central-difference schemes to upwind differences in areas, where the standard central differences do not provide monotone approximations for the convective terms, is used either explicitly or implicitly. Such regularized difference schemes for convection-diffusion equations with nondivergent and divergent convective terms are discussed in detail in the book [23] containing many references to various authors. If we consider 1D problems, then the simplest way to obtain unconditionally monotone schemes is based on a transformation of the governing equation, where the convective and diffusive transport are rewritten as the diffusive transport of an auxiliary variable. In this case, we arrive at exponential schemes proposed in the works [3, 24], which are in common use with many modifications in computational practice. Stability of difference schemes for time-dependent problems is investigated in Hilbert spaces of grid functions on the basis of the general theory of stability (correctness) for operator-difference schemes, which was developed by A. Samarskii [19,20]. The main results on stability of schemes for convection-diffusion problems can be found in [21, 23]. Different classes of unconditionally stable schemes of the first- and second-order accuracy in time are constructed in the most simple manner for the problems, where the convective terms are taken in the skew-symmetric form, i.e., as the average value of the nondivergent and divergent forms. For convection-diffusion problems with the convective terms in the nondivergent or divergent forms, it is natural to focus on the use of the Banach spaces L∞ (Ω) and L1 (Ω), respectively. Investigation of stability in Banach spaces L∞ (ω) and L1 (ω) (the discrete analogs of L∞ (Ω) and L1 (Ω)) is usually carried out on the basis of the maximum principle. In our work, we exploit the concept of the logarithmic norm. Using this mathematical tool, it is easy to establish stability properties of difference schemes for 1D convection-diffusion problems with central and upwind differences. Here monotonicity properties are studied for exponential schemes designed for convection-diffusion problems with the convective terms written both in the nondivergent and divergent forms. Here we consider the properties of stability and monotonicity for the schemes, which are well known in computational practice for solving the convection-diffusion equation. These schemes are not used in up-to-date CFD computations, but they provide the methodological basis for constructing new advanced schemes. Using the logarithmic norm technique, it is possible to consider in a simple and compact way two of the most important formulations of the convection-diffusion equation, where the convective terms are written in the divergent and nondivergent forms. Moreover, in this approach, there is no restriction on the compressibility of a medium. The properties of stability and monotonicity are studied on the basis of a unified methodology via selection of the corresponding Banach spaces of grid functions L∞ (ω) or L1 (ω). Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

Schemes for Unsteady Convection-Diffusion Problems

187

As it was mentioned above, the present study is restricted to the investigation of stability and monotonicity only. Other essential issues of constructing computational algorithms for the convection-diffusion equation are available elsewhere. For example, convergence of the approximate solution to the exact solution of a differential problem is discussed in the works [6,9]. A very important practical case of convection-dominated problems is studied in [17,25]. In the framework of the subjects of our theoretical study (sufficient conditions for stability and monotonicity), any numerical experiments are unnecessary. Applying coordinate-wise transformations, stable and monotone schemes are developed for multidimensional problems. Issues of constructing schemes of splitting with respect to spatial directions [12, 22] are studied separately. Stability conditions for monotone schemes are obtained in L∞ (ω) and L1 (ω). For more general problems of convection-diffusionreaction, the solution norm can grow exponentially with time. Possibilities of constructing special unconditionally stable schemes that take into account this feature are investigated, too.

1. Convection-Diffusion Equations Let us consider the class of model time-dependent convection-diffusion problems with a constant (independent of time, but depending on spatial coordinates) coefficient of diffusive transport. Coefficients of convective transport for applied problems depend on both space and time. To simplify the material presented here, we start with the 1D convection-diffusion problems. Consider the time-dependent convection-diffusion equation with convective terms in the nondivergent form ∂u ∂u ∂ ∂u + v(x, t) − k(x) = f (x, t) (1) ∂t ∂x ∂x ∂x for 0 < x < l,

t > 0.

Here v(x, t), k(x) and f (x, t) are given, sufficiently smooth, and bounded functions and k(x, t) > κ > 0. This equation is supplemented with the simplest homogeneous Dirichlet boundary conditions: u(0, t) = 0, u(l, t) = 0, t > 0. (2) In addition, the initial condition is specified as u(x, 0) = u0 (x),

0 < x < l.

(3)

The second important example is the unsteady equation of convection-diffusion with convective terms written in the divergent form ∂ ∂ ∂u ∂u + (v(x, t)u) − k(x) = f (x, t). (4) ∂t ∂x ∂x ∂x Consider the set of functions u(x, t) satisfying the boundary conditions (2). The transient problem of convection-diffusion may be written in the form of the operator-differential equation du + Au = f (t), A = A(t) = C(t) + D, (5) dt Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

188


where C(t) denotes the convective transport operator, and D stands for the operator of diffusive transport. The Cauchy problem for the evolutionary equation (5) is supplemented with the initial condition u(0) = u0 . (6) We recall some elementary a priori estimates for the convection-diffusion problems (1)– (3) and (2)–(4), which are derived from the maximum principle. The corresponding a priori estimates are obtained in the spaces L∞ (0, l) and L1 (0, l), where the norms are, respectively, Z l kvk∞ = max |v(x)|, kvk1 = |v(x)|dx. 0

m X

|aij |,

i = 1, 2, . . . , m

(14)

i = 1, 2, . . . , m

(15)

i6=j=1

(weak diagonal dominance by rows) or aii >

m X

|aji |,

i6=j=1

(weak diagonal dominance by columns). The logarithmic norm of a matrix A is defined [4, 7] by the number kE + δAk − 1 . δ→0+ δ

µ[A] = lim

For the logarithmic norm of a matrix in L∞ (consistent with (12)) and in L1 (consistent with (13)), we have the expressions m X µ∞ [A] = max aii + |aij | , 16i6m

m X µ1 [A] = max ajj + |aij | .

i6=j=1

16j6m

j6=i=1

In view of the restrictions (14) and (15), we have that the logarithmic norm of the matrix −A in the Cauchy problem (10)–(11) satisfies the inequality µ[−A] 6 0

(16)

in the corresponding space (in L∞ for (14), and in L1 for (15)). Among the properties of the logarithmic norm (see [4, 5]), we recall the following ones: 1. µ[cA] = cµ[A],

c = const > 0;

2. µ[cE + A] = c + µ[A],

c = const;

3. kAwk > max{−µ[−A], −µ[A]}kwk. Emphasis is placed on property 3, which allows to get easily the lower bound of the norm Aw. This bound can be combined with the standard upper bound of Aw: kAwk 6 kAkkwk. Let us study stability of difference schemes for solving the problem (10)–(11) numerically. We denote the approximate solution at the time level tn = nτ (where τ is a time step) as y n , and write the two-level difference scheme with weights y n+1 − y n + A(σy n+1 + (1 − σ)y n ) = ϕn , τ Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

(17)

190


where, e.g., A = A(σtn+1 + (1 − σ)tn ), with the initial data y 0 = u0 .

(18)

A sufficient condition for stability of the scheme (17)–(18) is formulated as the following statement. Theorem 2.1. Assume that in the Cauchy problem (10)–(11) the matrix A satisfies the restriction (14) (or (15)). Then the difference scheme with weights (17)–(18) is unconditionally stable for σ = 1, and it is conditionally stable for 0 6 σ < 1 in L∞ (in L1 ) if τ6

−1 1 max aii . 1 − σ 16i6m

(19)

Besides, the difference solution satisfies the a priori estimate ky n+1 k 6 ku0 k +

n X

τ kϕk k.

(20)

k=0

Proof. From (17), it follows that (E + στ A)y n+1 = (E − (1 − σ)τ A)y n + τ ϕn and therefore k(E + στ A)y n+1 k 6 k(E − (1 − σ)τ A)y n k + τ kϕn k.

(21)

For the left-hand side of (21), by the above-mentioned properties of the logarithmic norm and in view of (16), we have k(E + στ A)y n+1 k > −µ[−E − στ A] ky n+1 k = (1 − στ µ[−A])ky n+1 k > ky n+1 k. For the first term in the right-hand side of (21), we obtain k(E − (1 − σ)τ A)y n k 6 kE − (1 − σ)τ Ak ky n k. We investigate this estimate in more detail for L∞ . The case L1 is studied in a similar manner. Considering (12) and taking into account the condition of diagonal dominance (14), we have m X kE − (1 − σ)τ Ak = max 1 − (1 − σ)τ aii + aij 16i6m

i6=j=1

6 max |1 − (1 − σ)τ aii | + (1 − σ)τ 16i6m

m X

|aij |

i6=j=1

6 max |1 − (1 − σ)τ aii | + (1 − σ)τ aii 6 1 16i6m

with 0 6 σ 6 1 and under the restriction (19) on a time step. The substitution into (21) yields the inequality ky n+1 k 6 ky n k + τ kϕn k, which immediately implies the desired estimate (20) for stability with respect to the righthand side and the initial data. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


191

The above estimates for stability (20) in L∞ and L1 are directly associated with the monotonicity of the difference solution of the problem (17)–(18) under the assumption that the off-diagonal elements of the matrix A are non-positive. Let us prove the following statement. Theorem 2.2. Assume that in the schemes (17)–(18) the conditions of diagonal dominance (14) (or (15)) are fulfilled for aij 6 0,

i 6= j,

i, j = 1, 2, . . . , m

(22)

and let u0 > 0,

ϕn > 0,

n = 0, 1, . . . .

Then y n+1 > 0,

n = 1, 2, . . . ,

for any τ > 0 if σ = 1. But if 0 6 σ < 1, then this is true under the constraints on a time step (19). Proof. For the transition from the current time level to the next one, we have y n+1 + στ Ay n+1 = g n ,

n = 0, 1, . . . ,

(23)

where g n = y n − (1 − σ)τ Ay n + τ ϕn .

(24)

Suppose that y n > 0 (for n = 0 this is true from the assumptions of the theorem). We show that from this we have also the non-negativity of y n+1 (y n+1 > 0). We prove that under the assumptions of the diagonal dominance (14) (or (15)) and under the restrictions on a time step (19), for a non-negative y n and ϕn , we get g n > 0. In view of (24), we obtain gin

= (1 − (1 −

σ)τ aii )yin

− (1 − σ)τ

m X

aij yjn > (1 − (1 − σ)τ aii )yin > 0.

j6=i,j=1

In the conditions of the theorem, the matrix of the system of linear algebraic equations (23) is an M-matrix, i.e., we have: strong diagonal dominance, positive diagonal elements, and non-positive off-diagonal elements of the matrix. Because of this, from g n > 0, it follows that y n+1 > 0. Apply the derived results to studying stability and monotonicity of difference schemes for time-dependent problems of convection-diffusion in the nondivergent and divergent forms.

3. Difference Schemes for Convection-Diffusion Equations For simplicity, we restrict ourselves to uniform grids. Within the interval [0, l], we introduce a grid ω ¯ ≡ ω ∪ ∂ω = {x | x = xi = ih, i = 0, 1, . . . , N, N h = l}, where ω is the set of interior nodes: ω = {x | x = xi = ih, i = 1, 2, . . . , N − 1, N h = l}. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

192


After discretization in space of the model convection-diffusion problems with the homogeneous boundary conditions (1)–(3) and (2)–(4), we arrive at the problem (10)–(11), where m = N − 1 and the difference solution is defined as wi (t) = w(x, t), x ∈ ω. The discrete diffusion operator is specified, e.g., as follows: Dw = −

1 k(x + 0.5h)(w(x + h, t) − w(x, t)) h2 1 + 2 k(x − 0.5h)(w(x, t) − w(x − h, t)), h

x∈ω

(25)

with w(x, t) = 0,

x ∈ ∂ω.

(26)

Discretization of convective transport terms is conducted in such a way that v(x, t) are defined at the half-integer grid points ω ¯ . For operators of convective transport in the nondivergent form (equation (1)), in view of (26), we put 1 v(x + 0.5h, t)(w(x + h, t) − w(x, t)) 2h 1 + v(x − 0.5h, t)(w(x, t) − w(x − h, t)), x ∈ ω. (27) 2h A similar approximation of second order with respect to h for the convective transport operator in the divergent form (equation (4)) leads to Cw =

1 v(x + 0.5h, t)(w(x + h, t) + w(x, t)) 2h 1 (28) − v(x − 0.5h, t)(w(x, t) + w(x − h, t)), x ∈ ω. 2h Let us formulate a condition for stability and monotonicity of the schemes with weights (17)–(18) attributed to the problem (10)–(11), where Cw =

A=C +D

(29)

and D, C are specified according to (25)–(27) or (25)–(26), (28). Theorem 3.1. If σ = 1, then the difference scheme (17)–(18) with the relations (25)–(27), (29) (or with (25)–(26), (28)–(29)) is monotone, and the difference solution satisfies the a priori estimate (20) in L∞ (or in L1 ) under the restriction h|v(x ± 0.5h, t)| 6 2, k(x ± 0.5h)

x ∈ ω,

(30)

for any τ > 0. But if 0 6 σ < 1, then this is true under the constraint on a time step τ6

1 (1 − σ)γ

with 1 γ = max 2 k(x + 0.5h) + k(x − 0.5h) − x∈ω h for (27) and with 1 γ = max 2 k(x + 0.5h) + k(x − 0.5h) + x∈ω h for (28).

1 v(x + 0.5h, t) − v(x − 0.5h, t) 2h 1 v(x + 0.5h, t) − v(x − 0.5h, t) 2h

Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

(31)


193

Proof. Consider the case of the convection-diffusion equation (1)–(3) (approximations (25)– (27), (29)) in detail. The problem (2)–(4) (approximations (25)–(26), (28)–(29)) is examined in a similar way. To apply Theorems 2.1 and 2.2, write explicitly the elements of A. For (25)–(27), (29), we have 1 1 1 aii = 2 (ki+1/2 + ki−1/2 ) − vi+1/2 + vi−1/2 , h 2h 2h 1 1 ai,i−1 = − 2 ki−1/2 − vi−1/2 , h 2h 1 1 ai,i+1 = − 2 ki+1/2 + vi+1/2 , h 2h where ki±1/2 = k(x ± 0.5h), x ∈ ω. The condition of non-positivity of off-diagonal elements (22) holds for 1 1 k + vi−1/2 > 0, i−1/2 h2 2h

1 1 k − vi+1/2 > 0. i+1/2 h2 2h

(32)

In this case, diagonal dominance is assured. A spatial computational grid with the step from the conditions (30) satisfies the inequalities (32). Restrictions on a time step (19) are reduced to the particular condition (31). Thus, the conditions of Theorems 2.1 and 2.2 hold. This provides the stability and monotonicity of the difference solution of the convection-diffusion problem under the above restrictions on a time step. To overcomes restrictions (30) on a spatial grid, we apply upwind approximations for convective terms. Introduce notation v(x, t) = v + (x, t) + v − (x, t), 1 v + (x, t) = (v(x, t) + |v(x, t)|) > 0, 2 1 − v (x, t) = (v(x, t) − |v(x, t)|) 6 0. 2 Instead of (27), we put Cw =

1 − v (x + 0.5h, t)(w(x + h, t) − w(x, t)) h 1 + v + (x − 0.5h, t)(w(x, t) − w(x − h, t)). h

(33)

For the convective transport in the divergent form, we get Cw =

1 − v (x + 0.5h, t)w(x + h, t) − v − (x − 0.5h, t)w(x, t) h 1 + v + (x + 0.5h, t)w(x, t) − v + (x − 0.5h, t)w(x − h, t) . h

(34)

Theorem 3.2. If σ = 1, then the difference scheme (17)–(18) with the relations (25), (27), (29), (33) (or with (25)–(26), (29), (34)) is monotone, and the difference solution satisfies the a priori estimate (20) in L∞ (or in L1 ) for any τ > 0. But if 0 6 σ < 1, then this is true under the constraints on a time step (31) with 1 1 − + γ = max 2 k(x + 0.5h) + k(x − 0.5h) − v (x + 0.5h, t) − v (x − 0.5h, t) x∈ω h h Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

194


for (33) and with 1 1 γ = max 2 k(x + 0.5h) + k(x − 0.5h) + v + (x + 0.5h, t) − v − (x − 0.5h, t) x∈ω h h for (34). In particular, the fully implicit scheme (σ = 1) is unconditionally stable and monotone. The principal shortcomings of the above schemes are connected with upwind approximations for the convective terms (33) and (34) – these schemes demonstrate the first-order approximation in space. Schemes on the basis of the central-difference approximations (27) and (28) are more accurate – they have the second-order spatial approximation.

4. Exponential Schemes It is convenient to construct monotone schemes by means of transforming the original convection-diffusion equation, i.e., by eliminating the convective terms. The equation (1) may be written as 1 ∂ ∂u ∂u − k(x)χ(x, t) = f (x, t), (35) ∂t χ(x, t) ∂x ∂x where

Z χ(x, t) = exp −

x

0

v(s, t) ds . k(s)

(36)

The equation (4) is reduced to ∂u ∂ k(x) ∂(χ(x, t)u) − = f (x, t). ∂t ∂x χ(x, t) ∂x

(37)

Further, we can introduce discretizations in space. Similarly to (25), for the grid functions satisfying (26), it is possible to put in equation (35): Aw = −

1 h2 χ(x, t) +

k(x + 0.5h)χ(x + 0.5h, t)(w(x + h, t) − w(x, t))

1 h2 χ(x, t)

k(x − 0.5h)χ(x − 0.5h, t)(w(x, t) − w(x − h, t)),

where Z χ(x − 0.5h, t) = exp −

x−0.5h

0

v(s, t) ds . k(s)

Taking into account that Z χ(x − 0.5h, t) = χ(x, t) exp − x

x−0.5h

v(s, t) ds , k(s)

with a precision of O(h2 ), we put χ(x − 0.5h, t) = χ(x, t) exp(θ(x, t)h) Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

(38)


195

with notation

v(x, t) . 2k(x) Therefore, instead of (38), we can use the following approximation: θ(x, t) =

1 k(x + 0.5h) exp(−θ(x, t)h)(w(x + h, t) − w(x, t)) h2 1 + 2 k(x − 0.5h) exp(θ(x, t)h)(w(x, t) − w(x − h, t)). h For equation (37), similarly to (38), we put Aw = −

(39)

k(x + 0.5h) χ(x + h, t)w(x + h, t) − χ(x, t)w(x, t) + 0.5h, t) k(x − 0.5h) + 2 χ(x, t)w(x, t) − χ(x − h, t)w(x − h, t) . h χ(x − 0.5h, t)

Aw = −

h2 χ(x

Simplifying this expression, we obtain 1 k(x + 0.5h) exp(−θ(x + h, t)h)w(x + h, t) h2 1 + 2 k(x + 0.5h) exp(θ(x, t)h)w(x, t) h 1 + 2 k(x − 0.5h) exp(−θ(x, t)h)w(x, t) h 1 − 2 k(x − 0.5h) exp(θ(x − h, t)h)w(x − h, t). (40) h Using the above-introduced approximations for the convection-diffusion operator, we can construct monotone schemes. The primary statement is formulated as follows. Aw = −

Theorem 4.1. Assume that on the set of grid functions (26) the operator A is defined according to (39) (or (40)). In this case, if σ = 1, then the difference scheme (17)–(18) is monotone, and the difference solution satisfies the a priori estimate (20) in L∞ (or in L1 ) for any τ > 0. But if 0 6 σ < 1, then this is true under the constraints on a time step (31) with 1 γ = max 2 k(x + 0.5h) exp(−θ(x, t)h) + k(x − 0.5h) exp(θ(x, t)h) x∈ω h for (39) and 1 γ = max 2 k(x + 0.5h) exp(θ(x, t)h) + k(x − 0.5h) exp(−θ(x, t)h) x∈ω h for (40). Proof. In the case of (39), for the matrix elements, we have 1 k exp(−θ h) + k exp(θ h) , i i i+1/2 i−1/2 h2 1 ai,i−1 = − 2 ki−1/2 exp(θi h), h 1 ai,i+1 = − 2 ki+1/2 exp(−θi h). h Checking diagonal dominance by rows and the non-positivity of the off-diagonal elements is evident. aii =


196


In the case (40), we obtain 1 k exp(θ h) + k exp(−θ h) , i i i+1/2 i−1/2 h2 1 = − 2 ki−1/2 exp(θi−1 h), h 1 = − 2 ki+1/2 exp(−θi+1 h). h

aii = ai,i−1 ai,i+1

In view of the non-positivity of the off-diagonal elements, the condition of diagonal dominance by columns (15) takes the form aii > −ai−1,i − ai+1,i and it is obviously true. Thus, the conditions for stability and monotonicity are the same as for schemes with the upwind approximations of convective terms (Theorem 3.2). However, discretization in space is of second order as for schemes with the central-difference approximations (Theorem 3.1). Some complications in evaluating coefficients of the difference operator lead to a slight increasing of computational costs.

5. Multidimensional Problems Possibilities of constructing second-order monotone schemes for time-dependent equations of convection-diffusion are examined on model 2D problems. In a rectangle Ω = {x | x = (x1 , x2 ), 0 < xα < lα , α = 1, 2}, we study the unsteady convection-diffusion equation with the convective terms written in the nondivergent form 2

2

X ∂ ∂u ∂u ∂u X + vα (x, t) k(x) = f (x, t). − ∂t α=1 ∂xα α=1 ∂xα ∂xα

(41)

This equation is supplemented with homogeneous Dirichlet boundary conditions: u(x, t) = 0,

x ∈ ∂Ω,

t > 0.

(42)

In addition, it is supplemented with the initial condition u(x, 0) = u0 (x),

x ∈ Ω.

(43)

Similarly to the 1D case, the second example is the time-dependent equation of convection-diffusion with the convective transport written in the divergence form 2 2 X ∂u X ∂ ∂ ∂u + (vα (x, t)u) − k(x) = f (x, t). ∂t α=1 ∂xα ∂x ∂x α α α=1

(44)

The convection-diffusion operators in multidimensional problems are represented as the sum of 1D convection-diffusion operators. Because of this, in constructing monotone schemes Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


197

for multidimensional problems, we can apply the above approximations designed for the 1D operators of convection-diffusion. Similarly to (35) and (36), rewrite equation (41) as 2

1 ∂ ∂u ∂u X − k(x)χα (x, t) = f (x, t), ∂t α=1 χα (x, t) ∂xα ∂xα

(45)

where now Z χ1 (x, t) = exp − 0

x1

v1 (s, x2 , t) ds , k(s, x2 )

Z χ2 (x, t) = exp − 0

x2

v2 (x1 , s, t) ds . k(x1 , s)

(46)

A similar transformation for (44) yields 2

∂u X ∂ k(x) ∂(χα (x, t)u) − = f (x, t). ∂t α=1 ∂xα χα (x, t) ∂xα

(47)

For simplicity, we use a uniform grid in each spatial direction. For grids in separate directions xα , α = 1, 2, we use the notation introduced above: ω ¯ α = {xα | xα = iα hα , iα = 0, 1, . . . , Nα , Nα hα = lα }, and ωα = {xα | xα = iα hα , iα = 1, 2, . . . , Nα − 1, Nα hα = lα }. For the grid in Ω, we put ω ¯ ≡ ω ∪ ∂ω = ω ¯1 × ω ¯2,

ω = ω1 × ω2 .

After discretization in space of the boundary value problems (42)–(43), (45) and (42)– (43), (46), we arrive at the problem (10)–(11), where A = A1 + A2 ,

(48)

and Aα , α = 1, 2 are the 1D grid operators of convection-diffusion. On the set of grid functions such that w(x, t) = 0, x ∈ ∂ω (49) for equation (45), similarly to (39), we put A1 w = −

1 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x, t)h1 )w(x1 + h1 , x2 , t) h21 1 + 2 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x, t)h1 )w(x, t) h1 1 + 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x, t)h1 )w(x, t) h1 1 − 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x, t)h1 )w(x1 − h1 , x2 , t), h1 Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

(50)

198


A2 w = −

1 k(x1 , x2 + 0.5h2 ) exp(−θ2 (x, t)h2 )w(x1 , x2 + h2 , t) h22 1 + 2 k(x1 , x2 + 0.5h2 ) exp(−θ2 (x, t)h2 )w(x, t) h2 1 + 2 k(x1 , x2 − 0.5h2 ) exp(θ2 (x, t)h2 )w(x, t) h2 1 − 2 k(x1 , x2 − 0.5h2 ) exp(θ2 (x, t)h2 )w(x1 , x2 − h2 , t), h2

(51)

where

v1 (x, t) v2 (x, t) , θ2 (x, t) = , x ∈ ω. 2k(x) 2k(x) In the case of (47), we have (see (40)) 1 A1 w = − 2 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x1 + h1 , x2 , t)h1 )w(x1 + h1 , x2 , t) h1 1 + 2 k(x1 + 0.5h1 , x2 ) exp(θ1 (x, t)h)w(x, t) h1 1 + 2 k(x1 − 0.5h1 , x2 ) exp(−θ1 (x, t)h1 )w(x, t) h1 1 − 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x1 − h1 , x2 , t)h)w(x1 − h1 , x2 , t), h1 1 A2 w = − 2 k(x1 , x2 + 0.5h2 ) exp(−θ2 (x1 , x2 + h2 , t)h1 )w(x1 , x2 + h2 , t) h2 1 + 2 k(x1 , x2 + 0.5h2 ) exp(θ2 (x, t)h)w(x, t) h2 1 + 2 k(x1 , x2 − 0.5h2 ) exp(−θ2 (x, t)h1 )w(x, t) h2 1 − 2 k(x1 , x2 − 0.5h2 ) exp(θ2 (x1 , x2 − h2 , t)h)w(x1 , x2 − h2 , t). h2 Similarly to Theorem 4.1, the following statement is proved. θ1 (x, t) =

(52)

(53)

Theorem 5.1. Assume that on the set of grid functions (49) the operator A is defined according to (48), (50)–(51) (or (48), (52)–(53)). In this case, the difference scheme (17)– (18) is monotone, and the difference solution satisfies the a priori estimate (20) in L∞ (or in L1 ) for any τ > 0 if σ = 1. But if 0 6 σ < 1, then this is true under the constraints on a time step (31) with n1 1 γ = max 2 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x, t)h1 ) x∈ω h1 h1 o 1 1 + 2 k(x1 , x2 + 0.5h2 ) exp(−θ2 (x, t)h2 ) + 2 k(x1 , x2 − 0.5h2 ) exp(θ2 (x, t)h2 ) h2 h2 for (50)–(51) and with n1 1 γ = max 2 k(x1 + 0.5h1 , x2 ) exp(θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ) x∈ω h1 h1 o 1 1 + 2 k(x1 , x2 + 0.5h2 ) exp(θ2 (x, t)h2 ) + 2 k(x1 , x2 − 0.5h2 ) exp(−θ2 (x, t)h2 ) h2 h2 for (52)–(53). Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


199

6. Locally One-Dimensional Schemes Computational implementation of the exponential schemes (17)–(18), (48)–(51) and (17)– (18), (48)–(49), (52)–(53) involves the inversion of the non-selfadjoint elliptic grid operators E + στ A, where the matrix has strong diagonal dominance either by rows or by columns. To determine the numerical solution at a new time level, we can apply iterative methods. Another possibility is to use locally one-dimensional schemes, which are based on the splitting (48) [19,30]. Intending to 3D generalizations, we restrict ourselves to componentwise splitting schemes [12, 22]. Rewrite the difference equation (17) as follows: y n+1 = Sy n + τ ϕn ,

(54)

where S is the transition operator. For the scheme with weights (17), we have S = (E + στ A)−1 (E + (σ − 1)τ A).

(55)

For the scheme (18), (54) we have the stability condition kSk 6 1.

(56)

Monotonicity is ensured by the fact that the matrices (E + στ A)−1 and E + (σ − 1)τ A are M-matrices. Splitting schemes are constructed using transition operators for the individual terms in the additive representation (48). Let us define Sα (τ ) = (E + στ Aα )−1 (E + (σ − 1)τ Aα ),

α = 1, 2.

(57)

Instead of (55), we employ S = S1 (τ )S2 (τ ).

(58)

The stability condition (56) is true if kSα k 6 1,

α = 1, 2.

(59)

For the monotonicity of the scheme (54), (58), it is sufficient to require that the individual matrices Sα , α = 1, 2 will be M-matrices. For any value of σ, only the first-order accuracy with respect to τ is possible. Numerical implementation of the scheme (18), (54), (57)–(58) can be arranged using locally one-dimensional schemes with weights, i.e., y n+α/2 − y n+(α−1)/2 + Aα σy n+α/2 + (1 − σ)y n+(α−1)/2 = ϕnα , τ

α = 1, 2,

(60)

where, e.g., ϕn1 = 0,

ϕn2 = ϕn .

Theorem 6.1. Assume that on the set of grid functions (49) the operators Aα , α = 1, 2 are defined according to (50)–(51) (or (52)–(53)). In this case, if σ = 1, then the locally onedimensional difference scheme (18), (60) is monotone, and the difference solution satisfies Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

200


the a priori estimate (20) in L∞ (or in L1 ) for any τ > 0. But if 0 6 σ < 1, then this is true under the constraints on a time step (31) with n1 1 γ = max 2 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x, t)h1 ), x∈ω h1 h1 o 1 1 k(x , x + 0.5h ) exp(−θ (x, t)h ) + k(x , x − 0.5h ) exp(θ (x, t)h ) 1 2 2 2 2 1 2 2 2 2 h22 h22 for (50)–(51) or with n1 1 γ = max 2 k(x1 + 0.5h1 , x2 ) exp(θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ), x∈ω h1 h1 o 1 1 k(x , x + 0.5h ) exp(θ (x, t)h ) + k(x , x − 0.5h ) exp(−θ (x, t)h ) 1 2 2 2 2 1 2 2 2 2 h22 h22 for (52)–(53). Proof. Conditions for stability and monotonicity are verified for each individual equation of (60). In particular, when the operators Aα , α = 1, 2 are defined by (50)–(51), for the first equation, we have ky n+1/2 k 6 ky n k with γ = max x∈ω

o n1 1 k(x + 0.5h , x ) exp(−θ (x, t)h ) + k(x − 0.5h , x ) exp(θ (x, t)h ) . 1 1 2 1 1 1 1 2 1 1 h21 h21

For the second equation, we get ky n+1 k 6 ky n+1/2 k + τ kϕn k with γ = max x∈ω

n1 o 1 k(x , x + 0.5h ) exp(−θ (x, t)h ) + k(x , x − 0.5h ) exp(θ (x, t)h ) . 1 2 2 2 2 1 2 2 2 2 h22 h22

Monotonicity of locally one-dimensional schemes under consideration is established in a similar way. Other classes of splitting schemes can be applied, too. In this regard, we highlight the class of additively averaged schemes [22, 27]. Instead of a multiplicative representation of the transition operator (58), we can employ the additive representation 1 S = (S1 (2τ ) + S2 (2τ )) (61) 2 with preserving the first-order approximation in time for the scheme (54). For the scheme (54), (57), (61), we present another variant of numerical implementation. Define the auxiliary functions yαn+1 , α = 1, 2 from yαn+1 − yαn (62) + Aα (σyαn+1 + (1 − σ)yαn ) = 0. 2τ For the approximate solution at the new time level, we put 1 y n+1 = (y1n+1 + y2n+1 ) + τ ϕn . (63) 2 Conditions for stability and monotonicity for this additively averaged locally one-dimensional scheme are formulated in the following theorem. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


201

Theorem 6.2. Assume that on the set of grid functions (49) the operators Aα , α = 1, 2 are defined according to (50)–(51) (or (52)–(53)). In this case, if σ = 1, then the additively averaged locally one-dimensional scheme (18), (62)–(63) is monotone, and the difference solution satisfies the a priori estimate (20) in L∞ (or in L1 ) for any τ > 0. But if 0 6 σ < 1, then this is true under the constraints on a time step (31) with n1 1 γ = 2 max 2 k(x1 + 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(θ1 (x, t)h1 ), x∈ω h1 h1 o 1 1 k(x1 , x2 + 0.5h2 ) exp(−θ2 (x, t)h2 ) + 2 k(x1 , x2 − 0.5h2 ) exp(θ2 (x, t)h2 ) h22 h2 for (50)–(51) and with n1 1 γ = 2 max 2 k(x1 + 0.5h1 , x2 ) exp(θ1 (x, t)h1 ) + 2 k(x1 − 0.5h1 , x2 ) exp(−θ1 (x, t)h1 ), x∈ω h1 h1 o 1 1 k(x1 , x2 + 0.5h2 ) exp(θ2 (x, t)h2 ) + 2 k(x1 , x2 − 0.5h2 ) exp(−θ2 (x, t)h2 ) h22 h2 for (52)–(53). Additively average schemes, on the one hand, demonstrate lower accuracy in comparison with schemes of componentwise splitting, but on the other hand, they are more promising in terms of parallel computing – the components yαn+1 , α = 1, 2 are determined (see (62)) independently of each other.

7. Convection-Diffusion-Reaction Problems Among the most important and general problems, special attention should be given to problems of convection-diffusion-reaction. Additional low-order terms in the parabolic equation, which describe chemical reactions, result in the fact that the solution norm for the homogeneous equation can exponentially grow with time. Such a behavior of the solution leads to the necessity to consider %-stability of difference schemes [19–21]. Different nature of convective, diffusive, and reactive transport can be taken into account in constructing approximations in time. The most pronounced inhomogeneity of approximation in time is expressed in explicit-implicit schemes. In this case, for solving transient problems numerically, some part of terms of the problem operator is approximated by explicit relationships, whereas other part is discretized via implicit formulas. Explicitimplicit schemes are widely used in the numerical solution of convection-diffusion problems [1]. In view of subordination of the convective transport operator to the diffusion operator, it is possible to prove unconditional stability of explicit-implicit schemes for unsteady convection-diffusion problems [23]. Such techniques are also applied to study problems of diffusion-reaction. In this case (see, e.g., [18]), the diffusive transport is treated implicitly, whereas explicit approximations are used to account reactions (source terms). In [28], schemes for convection-diffusion-reaction problems are constructed, which are unconditionally stable in Hilbert spaces of grid functions, with the operator splitting, where one part of the reaction operator is taken from a lower time level, and the other part is referred to the next (upper) time level. Here we use a distinct approach to the construction of unconditionally stable schemes for solving unsteady convection-diffusion-reaction equations numerically, which in some sense is Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

202


similar to nonstandard schemes. To take into account peculiarities of the transient problems under consideration, we modify the standard two-level schemes with weights. Some examples of such nonstandard discretizations in time can be found, e.g., in [13, 14]. We restrict ourselves to simple 1D problems, where, similarly to (1) and (4), we seek the solution of the boundary value problem for the convection-diffusion-reaction equation with convective transport written in the nondivergent form ∂u ∂ ∂u ∂u + r(x, t)u + v(x, t) − k(x) = f (x, t), (64) ∂t ∂x ∂x ∂x or the convective terms in the divergent form ∂u ∂ ∂ ∂u + r(x, t)u + (v(x, t)u) − k(x) = f (x, t). ∂t ∂x ∂x ∂x

(65)

The reaction coefficient r(x, t) is considered under the constraint r(x, t) > µ,

0 < x < l,

µ = const < 0.

(66)

The simplest way to derive an a priori estimate for the boundary value problems (2)–(3), (64), (66) and (2)–(3), (65)–(66) is based on introducing a new unknown variable. We set u(x, t) = exp(−µt)z(x, t).

(67)

In this case, for z(x, t), from (2)–(3), (64) and (67), we obtain the problem ∂z ∂w ∂ ∂z + (r(x, t) − µ)z + v(x, t) − k(x) = exp(µt)f (x, t), ∂t ∂x ∂x ∂x z(0, t) = 0, z(l, t) = 0, t > 0, z(x, 0) = u0 (x), 0 < x < l. The solution of this problem satisfies the estimate of type (7): Z t 0 exp(µθ)kf (x, θ)k∞ dθ. kz(x, t)k∞ 6 ku (x)k∞ + 0

In view of (67), we arrive at the estimate 0

Z

ku(x, t)k∞ 6 exp(−µt)ku (x)k∞ +

t

exp(−µ(t − θ))kf (x, θ)k∞ dθ

(68)

0

for the solution of the problem (2)–(3), (64), (66). Similarly, for (2)–(3), (65)–(66), we have Z t 0 ku(x, t)k1 6 exp(−µt)ku (x)k1 + exp(−µ(t − θ))kf (x, θ)k1 dθ. (69) 0

Discretization in space yields the Cauchy problem (10)–(11), where A = C + D + R.

(70)

Here the operators D, C are defined according to (25)–(27) or (25)–(26), (28), and for R, we set Rw = r(x, t)w, x ∈ ω. (71) To construct discretization in time, we use a transformation of the problem. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


203

Similarly to (67), in the problem (10)–(11), (66), (70), we put w(x, t) = exp(−µt)z(x, t),

x ∈ ω.

For z(x, t), we obtain the problem dz ˜ + A(t)z = exp(µt)φ(t), dt z(0) = u0 .

(72) (73)

The operator A˜ can be represented as A˜ = C + D + R − µE,

(74)

where E is the identity operator. Additional terms of the operator R − µE in (74) have (see (66), (71)) diagonal dominance both by rows and by columns. Similarly to (17)–(18), we solve the problem (72)–(73) using the two-level scheme with weights sn+1 − sn ˜ n+1 + (1 − σ)sn ) = ϕ˜n , + A(σs τ s0 = u0 ,

(75) (76)

where ˜ n+1 + (1 − σ)tn ), A˜ = A(σt

ϕ˜n = exp(µστ ) exp(µtn )ϕ(σtn+1 + (1 − σ)tn ).

Conditions for stability of the scheme (75)–(76) are checked similarly to the stability conditions of (17)–(18). The estimate for stability with respect to the initial data and the right-side has the form n X n+1 0 τ kϕ˜k k. (77) ks k 6 ku k + k=0

To construct the scheme with weights for (10)–(11), (66), (70), we apply relations sn+1 = exp(µtn+1 )y n+1 ,

sn = exp(µtn )y n .

Substitution into (75)–(76) yields the scheme exp(µτ )y n+1 − y n + A˜ σ exp(µτ )y n+1 + (1 − σ)y n = exp(µστ )ϕn , τ y 0 = u0 .

(78) (79)

From (77), it follows that the estimate ky

n+1

n+1

k 6 exp(−µt

0

)ku k +

n X

τ exp(µστ ) exp(−µ(tn+1 − tk ))kϕk k

k=0

holds. Thus, we obtain the following result. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM

(80)

204


Theorem 7.1. If σ = 1, then the difference scheme (78)–(79) with relations (25)–(27), (70)– (71) (or with (25)–(26), (28), (70)–(71)) is monotone, and the difference solution satisfies the a priori estimate (80) in L∞ (or in L1 ) under the restriction (30) for any τ > 0. But if 0 6 σ < 1, then this is true under the restriction on a time step (31) with 1 1 v(x + 0.5h, t) − v(x − 0.5h, t) γ = max 2 k(x + 0.5h) + k(x − 0.5h) − x∈ω h 2h for (27) and with 1 1 γ = max 2 k(x + 0.5h) + k(x − 0.5h) + v(x + 0.5h, t) − v(x − 0.5h, t) x∈ω h 2h for (28).

References [1] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton, Implicit-explicit methods for time-dependent partial differential equations, SIAM J. Numer. Anal., 32 (1995), pp. 797–823. [2] G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 2000. [3] D. N. de G. Allen and R. V. Southwell, Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder, Quart. J. Mech. Appl. Math., 8 (1955), pp. 129–145. [4] K. Dekker and J. G. Verwer, Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984. [5] C. Desoer and H. Haneda, The measure of a matrix as a tool to analyze computer algorithms for circuit analysis, IEEE Trans. Circuit Theory, 19 (1972), pp. 480–486. [6] C. Grossmann, H.-G. Roos, and M. Stynes, Numerical Treatment of Partial Differential Equations, Springer, Berlin, 2007. [7] E. Hairer, S. P. Norsett, and G. Wanner, Solving Ordinary Differential Equations. I. Nonstiff Problems, Springer, Berlin, 1987. [8] W. Hundsdorfer and J. G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer, Berlin, 2003. [9] P. Knabner and L. Angermann, Numerical Methods for Elliptic and Parabolic Partial Differential Equations, Springer, New York, 2003. [10] O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Ural’ceva, Linear and Quasi-Linear Equations of Parabolic Type, Transl. Math. Monogr., 23, American Mathematical Society, Providence, 1968. [11] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics, vol. 6: Fluid Mechanics, 2nd edn., Butterworth-Heinemann, Oxford, 1987. [12] G. I. Marchuk, Splitting and alternating direction methods, in: Handbook of Numerical Analysis, vol. I, North-Holland, Amsterdam, 1990, pp. 197–462. [13] R. E. Mickens, Nonstandard finite difference schemes for differential equations, J. Difference Equ. Appl., 8 (2002), pp. 823–847. [14] R. E. Mickens (ed.), Advances in the Applications of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2005. [15] K. W. Morton, Numerical Solution of Convection-Diffusion Problems, Chapman & Hall London, 1996. [16] M. H. Protter and H. F. Weinberger, Maximum Principles in Differential Equations, Springer, Berlin, 1984. Brought to you by | Sterling Memorial Library Authenticated | 142.58.101.27 Download Date | 10/6/13 4:02 AM


205

[17] H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations Convection-Diffusion-Reaction and Flow Problems, Springer, Berlin, 2008. [18] S. J. Ruuth, Implicit-explicit methods for reaction-diffusion problems in pattern formation, J. Math. Biol. 34 (1995), pp. 148–176. [19] A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001. [20] A. A. Samarskii and A. V. Gulin, Stability of Difference Schemes (in Russian), Nauka, Moscow, 1973. [21] A. A. Samarskii, P. P. Matus, and P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer, Dordrecht, 2002. [22] A. A. Samarskii and P. N. Vabishchevich, Additive Schemes for Problems of Mathematical Physics (in Russian), Nauka, Moscow, 1999. [23] A. A. Samarskii and P. N. Vabishchevich, Numerical Methods for Solving Convection-Diffusion Problems (in Russian), URSS, Moscow, 1999. [24] D. L. Scharfetter and H. K. Gummel, Large-signal analysis of a silicon read diode oscillator, IEEE Trans. Electron Devices, ED-16 (1969), pp. 4–77. [25] G. I. Shishkin and L. P. Shishkina, Difference Methods for Singular Perturbation Problems, CRC Press, Boca Raton, 2008. [26] J. C. Tannehill, D. A. Anderson, and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Taylor & Francis, Washington, DC, 1997. [27] P. N. Vabishchevich, Construction of splitting schemes based on transition operator approximation, Comput. Math. Math. Phys., 52 (2012), pp. 235–244. [28] P. N. Vabishchevich and M. V. Vasil’eva, Explicit-implicit schemes for convection-diffusion-reaction problems, Siberian J. Num. Math., 15 (2012), pp. 359–369. [29] P. Wesseling, Principles of Computational Fluid Dynamics, Springer, Berlin, 2001. [30] N. N. Yanenko, The Method of Fractional Steps. The Solution of Problems of Mathematical Physics in Several Variables, Springer, New York, 1971.
