PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 91(105) (2012), 155–162 DOI: 10.2298/PIM1205155V

DOMAIN DECOMPOSITION SCHEMES FOR THE STOKES EQUATION

Petr N. Vabishchevich

Abstract. Numerical algorithms for solving problems of mathematical physics on modern parallel computers employ various domain decomposition techniques. Domain decomposition schemes are developed here to solve numerically initial/boundary value problems for the Stokes system of equations in the primitive variables pressure–velocity. Unconditionally stable schemes of domain decomposition are based on the partition of unit for a computational domain and the corresponding Hilbert spaces of grid functions.

1. Introduction

In computational fluid dynamics [1, 10] there are employed numerical algorithms based on using the primitive variables pressure-velocity. The main difficulties in this approach are connected with the calculation of the pressure. In studying transient problems the corresponding elliptic Neumann problem for the pressure is derived as the result of employment of one or another scheme of splitting with respect to physical processes [3, 4].

Domain decomposition methods are used for the numerical solution of boundary value problems for partial differential equations on parallel computers. They are in most common use for stationary problems [7, 11]. Computational algorithms with and without overlapping of subdomains are employed in this case in synchronous (sequential) and asynchronous (parallel) algorithms.

For transient problems it seems to be more suitable to utilize iteration-free variants of domain decomposition techniques [6, 8] which are best suited to peculiarities of a problem (evolution in time). In these regionally-additive schemes a transition to a new time level is performed via solving problems in particular subdomains.

The regionally-additive schemes for the Navier–Stokes equations in the primitive variables are discussed in [2]. In simulation of incompressible flows an elliptic problem for the pressure can be changed to separate elliptic problems for the pressure in particular subdomains. Therefore, it is possible to construct iterationfree regionally-additive schemes for the Navier–Stokes equations. In this paper we propose a general approach to construct domain decomposition schemes for timedependent systems of equations. Using the partition of unit for a computational domain and the corresponding Hilbert spaces of grid functions we perform a transition to finding the individual components of the solution in the subdomains. The unsteady Stokes equations for an incompressible fluid is considered as a typical problem.

2010 Mathematics Subject Classification: Primary 65M06, 65M12; Secondary 76D07.

Key words and phrases: Viscous incompressible flows, numerical methods, domain decomposition techniques, operator-splitting schemes.

Dedicated to Academician Anton Bilimovic.

2. Stokes equations

Assume that the linear approximation is valid to describe a flow of an incompressible viscous fluid. In a region Ω with solid boundaries we can write equations of motion and continuity in the primitive variables pressure, velocity as follows

(2.1) \frac{∂u}{∂t} + gradp − νΔu = f(x, t),

divu = 0, x ∈ Ω, 0 < t ⩽ T. (2.2)

Here u is the velocity, p is the pressure, ν is the kinematic viscosity and Δ = divgrad is the Laplace operator. Equations (2.1), (2.2) are supplemented with the following condition for the single-valued evaluation of the pressure

(2.3)  ∫_{Ω} p(x, t) dx = 0, 0 < t ⩽ T.

No-slip, no-permeability conditions are specified on solid boundaries

(2.4)  u(x, t) = 0, x ∈ ∂Ω, 0 < t ⩽ T.

Some initial condition is also given

(2.5)  u(x, 0) = v(x), x ∈ Ω.

Let us rewrite problem (2.1)–(2.5) in an operator formulation. On the set of functions satisfying (2.3),2.4), we have the Cauchy problem

\frac{du}{dt} + Au + Bp = f, (2.6)

B^{∗}u = 0, 0 < t ⩽ T, (2.7)

u(0) = v. (2.8)

For these operators in the space L_{2}(Ω) we have A = A^{∗} ⩾ δE, δ = δ(Ω) > 0, where E is the unit (identity) operator. Adjointness of operators B = grad(B : L_{2}(Ω) → L_{2}(Ω)) and B^{∗} = − div(B^{∗} : L_{2}(Ω) → L_{2}(Ω)) follows from

∫_{Ω} u gradp dx + ∫_{Ω} divu p dx = 0.

For problem (2.6)–(2.8) the following simple a priori estimate is valid

(2.9)  ∥u(t)∥^{2} ⩽ ∥v∥^{2} + \frac{1}{2δ} ∫_{0}^{t} ∥f(θ)∥^{2}dθ,

where ∥ · ∥ is the norm in L_{2}(Ω). Estimate (2.9) will be for us a reference point when considering discrete problems.

3. Approximation in space

In this paper the main attention is paid to computational algorithms for the transition to a new time level, i.e., approximation in time. To construct discretization in time, operator-splitting schemes are used that allow to formulate a problem for the pressure in the most natural way. The problem of approximation in space is solved in the standard manner.

There are employed various types of grids: the non-staggered (collocated) grid, where both the pressure and velocity components are referred to the same points; next, partially staggered (ALE-type) grid, where the pressure is referred to an individual grid shifted in all space directions on a one-half step from the basic grid where all velocity components are defined; and finally, the staggered (MAC-type) grid, where the pressure is defined at the center points of grid cells, whereas the velocity components are referred to the corresponding faces of the cell.

For simplicity we consider here uniform rectangular non-staggered grids. Problem (2.1)–(2.5) is solved in a rectangle

Ω = {x | x = (x_{1}, x_{2}), 0 < x_{α} < l_{α} , α = 1, 2}.

The approximate solution is calculated at the points of a uniform rectangular grid in Ω:

ω̄ = {x | x = (x_{1}, x_{2}), x_{α} = i_{α} h_{α} , i_{α} = 0, 1, …, N_{α} , N_{α} h_{α} = l_{α} }

and let ω be the set of internal nodes (ω̄ = ω ∪ ∂ω).

For vector grid functions u(x) = 0, x ∈ ∂ω we define a Hilbert space H = L_{2}(ω) with the scalar product and norm

(u, v) = ∑_{x∈ω} u(x)v(x)h_{1}h_{2}, ∥u∥ = (u, u)^{1/2} .

The grid operator A is taken in the form A = −νΔ_{h}, where Δ_{h} is the grid Laplace operator:

Δ_{h}y = −\frac{1}{h_{1}^{2}} (y(x_{1} + h_{1}, x_{2}) − 2y(x_{1}, x_{2}) + y(x_{1} − h_{1}, x_{2}))

− \frac{1}{h_{2}^{2}} (y(x_{1}, x_{2} + h_{2}) − 2y(x_{1}, x_{2}) + y(x_{1}, x_{2} − h_{2})).

In H the operator A is selfadjoint and positive definite:

(3.1)  A = −νΔ_{h} = A^{∗} ⩾ νδ_{h}E, δ_{h} = ∑_{α=1}^{2} \frac{4}{h_{α}^{2}} sin^{2} \frac{πh_{α}}{2l_{α}} .

The pressure gradient is approximated by directed differences with an error O(h). We set B = grad_{h} at

(3.2)  Bp = {(Bp)_{1}, (Bp)_{2}}, x ∈ ω,

where

(Bp)_{1} = \frac{1}{h_{1}} (p(x_{1} + h_{1}, x_{2}) − p(x_{1}, x_{2})),

(Bp)_{2} = \frac{1}{h_{2}} (p(x_{1}, x_{2} + h_{2}) − p(x_{1}, x_{2})).

The set of points for the pressure evaluation is denoted as ω_{p} (ω_{p} ⊂ ω̄). For the grid divergence operator B^{∗} = − div_{h} we have

(3.3)  B^{∗}u = −\frac{1}{h_{1}} (u_{1}(x_{1}, x_{2}) − u_{1}(x_{1} − h_{1}, x_{2}))

− \frac{1}{h_{2}} (u_{2}(x_{1}, x_{2}) − u_{2}(x_{1}, x_{2} − h_{2})), x ∈ ω_{p}.

The adjointness property of the grid gradient and divergence operators is a consequence of the discrete equation

∑_{x∈ω} Bp(x)u(x)(x)h_{1}h_{2} + ∑_{x∈ω_{p}} p(x)B^{∗}u(x)(x)h_{1}h_{2} = 0,

which takes place on the set of vector grid functions u(x) = 0, x ∈ ∂ω.

In view of (3.1)–(3.3) after the spatial approximation of problem (2.6)–(2.8) we obtain the following problem

\frac{du}{dt} + Au + Bp = f, (3.4)

B^{∗}u = 0, 0 < t ⩽ T, (3.5)

u(0) = v. (3.6)

For the solution of problem (3.4)–(3.6) a priori estimate (2.9) holds, where now ∥·∥ is the norm in H = L_{2}(ω).

4. Domain decomposition

Let Ω be a combination of m particular subdomains Ω = Ω_{1} ∪ Ω_{2} ∪ ⋯∪ Ω_{m}. Particular subdomains can overlap one onto another. We shall construct the schemes of decomposition where the solution at the new time level for the initial problem is reduced to the sequential solution of problems in particular subdomains.

Let us define functions for domain Ω

(4.1)  η_{α} (x) = \{ \begin{matrix}> 0, & x ∈ Ω_{α} , \\ 0, & x ∉ Ω_{α} , \end{matrix} α = 1, 2, …, m.

Generally, see for example,[6, 8], domain decomposition schemes for unsteady problems are based on the partition of unit for the region Ω, where

∑_{α=1}^{m} η_{α} (x) = 1, x ∈ Ω.

It is more convenient to use a somewhat different partition where

(4.2)  ∑_{α=1}^{m} η_{α}^{2}(x) = 1, x ∈ Ω.

For the decomposition of computational domain (4.1), (4.2) we consider the following additive representation of the identity operator E in H = L_{2}(ω):

(4.3)  E = ∑_{α=1}^{m} χ_{α}^{2} , χ_{α} = η_{α} (x)E, x ∈ ω, α = 1, 2, …, m.

Taking into account (4.3) we have

(4.4)  u = ∑_{α=1}^{m} u_{α} , u_{α} = χ_{α} u, α = 1, 2, …, m.

To formulate an appropriate system of equations for determining components of the solution u_{α}, α = 1, 2, …, m, we multiply both sides of equation (3.4) by χ_{α}. This gives

(4.5)  \frac{du_{α}}{dt} + χ_{α} A ∑_{β=1}^{m} χ_{β}u_{β} + B_{α} p = f_{α} ,

where B_{α} = χ_{α} B, f_{α} = χ_{α} f, α = 1, 2, …, m. Taking into account that B_{α}^{∗} = B^{∗}χ_{α}, equation (3.5) in the new notation is written as

(4.6)  ∑_{α=1}^{m} B_{α}^{∗}u_{α} = 0, 0 < t ⩽ T.

The system of equations (4.5), (4.6) is supplemented by the initial conditions

(4.7)  u_{α} (0) = v_{α} , v_{α} = χ_{α} v, α = 1, 2, …, m.

For U = {u_{1}, u_{2}, …, u_{m}}, we define in the space H^{m} the norm and inner product as follows (U, V)_{m} = ∑ _{α=1}^{m} (u_{α} , v_{α}), ∥U∥_{m} = (U, U)_{m}^{1/2} . Let us multiply scalarly individual equations (4.5) by u_{α}, α = 1, 2, …, m and add them together. Next, multiply equation (4.6) by p. Taking into account (4.4), we obtain

\frac{1}{2} \frac{d}{dt} ∑_{α=1}^{m} (u_{α} , u_{α}) + (Au, u) = ∑_{α=1}^{m} (f_{α} , u_{α}).

This implies the a priori estimate

(4.8)  ∥U∥_{m}^{2} ⩽ exp(t)∥V ∥_{m}^{2} + ∫_{0}^{t} exp(t − θ)∥F (θ)∥_{m}^{2}dθ

for problem (4.5)–(4.7) with F = {f_{1}, f_{2}, …, f_{m}}.

5. Splitting scheme

In the construction of domain decomposition schemes we shall proceed from the scheme of splitting with respect to physical processes for the Cauchy problem (2.6)–(2.8). We shall use a simple additive scheme componentwise splitting [5, 9]. Let u^{n} be the difference solution at the time moment t^{n} = nτ, where τ = T/N > 0 is the time-step. Let us separate out a particular stage connected with the pressure impact[3, 4]. Thus, in the first stage we have

(5.1)  \frac{u^{n+1/2} − u^{n}}{τ} + Au^{n+1/2} = f^{n+1/2} .

The pressure gradient is treated only in the second stage:

\frac{u^{n+1} − u^{n+1/2}}{τ} + Bp^{n+1} = 0, (5.2)

B^{∗}u^{n+1} = 0. (5.3)

Implementation of (5.2), (5.3) consists of two steps. In the first step we solve the following problem for the pressure B^{∗}Bp^{n+1} = \frac{1}{τ}B^{∗}u^{n+1/2} , whereas in the second one we update the velocity: u^{n+1} = u^{n+1/2} − τBp^{n+1} . Multiplying (5.1) by u^{n+1/2} , we obtain

∥u^{n+1/2} ∥^{2} ⩽ ∥u^{n} ∥^{2} + \frac{τ}{δ_{h}} ∥f^{n+1/2} ∥^{2}.

Similarly, from (5.2) taking into account (5.3) we have ∥u^{n+1} ∥^{2} ⩽ ∥u^{n+1/2} ∥^{2}. Thus, we obtain the grid analog of estimate (2.9)

∥u^{n+1} ∥^{2} ⩽ ∥u^{n} ∥^{2} + \frac{τ}{2δ_{h}} ∥f^{n+1/2} ∥^{2}.

for difference scheme (5.1)–(5.3).

For simplicity, we shall construct splitting schemes for problem (4.5)–(4.7) by analogy with splitting scheme (5.1)–(5.3). The first half-step (viscous dissipation) is associated with the solution of equations

\frac{du_{α}}{dt} + χ_{α} A ∑_{β=1}^{m} χ_{β}u_{β} = f_{α} , α = 1, 2, …, m, t^{n} < t ⩽ t^{n+1/2} .

Taking into account the fact that the numerical solution is implemented via solving individual problems in the subdomains, the transition from time level t^{n} to level t^{n+1/2} can be realized as follows:

(5.4) \frac{u_{α}^{n+1/4} − u_{α}^{n}}{τ} + χ_{α} A ∑_{β=1}^{α−1} χ_{β}u_{β}^{n+1/4} + \frac{1}{2}χ_{α} Aχ_{α} u_{α}^{n+1/4} = f_{α}^{n+1/2} ,

α = 1, 2, …, m,

(5.5) \frac{u_{α}^{n+1/2} − u_{α}^{n+1/4}}{τ} + \frac{1}{2}χ_{α} Aχ_{α} u_{α}^{n+1/2} + χ_{α} A ∑_{β=α+1}^{m} χ_{β}u_{β}^{n+1/2} = 0,

α = 1, 2, …, m.

In view of (5.4), (5.5) in each subdomain Ω_{α} , α = 1, 2, …, m we must invert grid selfadjoint elliptic operator D_{α} = E + \frac{1}{2}χ_{α} Aχ_{α} for finding u_{α}^{n+1/4} (α = 1, 2, …, m) and u_{α}^{n+1/2} (α = m, m − 1, …, 1). In this case, outside subdomains Ω_{α} , α = 1, 2, …, m explicit calculations are used.

Stability of scheme (5.4), (5.5) will be investigated in H^{m}. Consider the operator

(5.6)  𝔸 = {A_{αβ} }, A_{αβ} = χ_{α} Aχ_{β}, α, β = 1, 2, …, m.

Taking into account (3.1), (4.2), we have 𝔸 = 𝔸^{∗} ⩾ 0 in H^{m}. Scheme (5.4), (5.5) is based on using the triangular splitting

(5.7)  𝔸 = 𝔸_{1} + 𝔸_{2}, 𝔸_{1} = 𝔸_{2}^{∗}.

Using notation (5.6), (5.7) we rewrite (5.4), (5.5) in the form

(5.8)  \frac{U^{n+1/4} − U^{n}}{τ} + 𝔸_{1}U^{n+1/4} = F^{n+1/2} ,
(5.9)  \frac{U^{n+1/2} − U^{n+1/4}}{τ} + 𝔸_{2}U^{n+1/2} = 0.

Taking into account that 𝔸_{α} ⩾ 0, α = 1, 2 in H^{m}, for (5.9) we immediately have

(5.10)  ∥U^{n+1/2} ∥_{m}^{2} ⩽ ∥U^{n+1/4} ∥_{m}^{2}.

Multiplying (5.8) by U^{n+1/4} , we obtain ∥U^{n+1/4} ∥_{m}^{2} ⩽ ∥U^{n}∥_{m}^{2}+2τ(F^{n+1/2} , U^{n+1/4} )_{m}. For the last term on the right-hand side we use the estimate

2τ(F^{n+1/2} , U^{n+1/4} )_{m} ⩽ (1 − exp(−τ))∥U^{n+1/4} ∥_{m}^{2} + \frac{τ^{2}}{1 − exp(−τ))} ∥F^{n+1/2} ∥_{m}^{2}.

This leads to the estimate

(5.11)  ∥U^{n+1/4} ∥_{m}^{2} ⩽ exp(τ)∥U^{n} ∥_{m}^{2} + τ∥F^{n+1/2} ∥_{m}^{2}.

The second half-step results from the pressure and is connected with the system of equations

\frac{du_{α}}{dt} + B_{α} p = 0, α = 1, 2, …, m,
∑_{α=1}^{m} B_{α}^{∗}u_{α} = 0, t^{n+1/2} < t ⩽ t^{n+1} .

Approximation in time for such systems were considered in [12]. We shall use the additive scheme

u_{α}^{n+1/2+β/2m} = u_{α}^{n+1/2+(β−1)/2m} , β ≠ α, β = 1, 2, …, m, (5.12)

\frac{u_{α}^{n+1/2+α/2m} − u_{α}^{n+1/2+(α−1)/2m}}{τ} + B_{α} p^{n+1/2+α/2m} = 0, (5.13)

B_{α}^{∗}u_{α}^{n+1/2+α/2m} = 0, α = 1, 2, …, m. (5.14)

The implementation of additive scheme (5.12)–(5.14) is conducted by analogy with scheme (5.2), (5.3).

For (5.12)–(5.14) we have ∥u_{α}^{n+1} ∥ ⩽ ∥u_{α}^{n+1/2} ∥, α = 1, 2, …, m and therefore

(5.15)  ∥U^{n+1} ∥_{m}^{2} ⩽ ∥U^{n+1/2} ∥_{m}^{2}.

Taking into account (5.10), (5.11) and (5.15) we obtain the desired stability estimate of the additive operator-difference scheme (5.4), (5.5), (5.12)–(5.14)

(5.16)  ∥U^{n+1} ∥_{m}^{2} ⩽ exp(τ)∥U^{n} ∥_{m}^{2} + τ∥F^{n+1/2} ∥_{m}^{2}.

Estimate (5.16) is the grid analog of (4.8) for the differential problem.

This allows us to formulate the following main result.

Theorem 5.1. The additive scheme of domain decomposition (5.4), (5.5), (5.12)–(5.14) is unconditionally stable and estimate (5.16) holds for the numerical solution.

References

1. J. D. Anderson and J. F. Wendt, Computational Fluid Dynamics, McGraw-Hill, 1995.

2. A. G. Churbanov and P. N. Vabishchevich, Domain decomposition techniques to solve the navier-stokes equations; in: B. Chetverushkin, A. Ecer, J. Periaux, N. Satofitka, and P. Fox (eds.), Parallel Computational Fluid Dynamics. Advanced Numerical Methods, Software and Applications, Amsterdam: Elsevier, 2004, pp. 145–152.

3. R. Glowinski, Finite element methods for incompressible viscous flow, in: Ciarlet, P. G. (ed.) et al., Numerical Methods For Fluids (Part 3), Amsterdam: North-Holland; Handb. Numer. Anal. 9, 2003.

4. J. L. Guermond, P. Minev, and J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (44–47) (2006), 6011–6045.

5. G. I. Marchuk, Splitting and alternating direction methods; in: Ciarlet, P. G. (ed.) et al., Handbook of Numerical Analysis 1, Amsterdam: North-Holland, 1990, 197-462.

6. T. P. A. Mathew, Domain Decomposition Methods for the Numerical Solution of Partial Differential Equations, Springer-Verlag, 2008.

7. A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford University Press, USA, 1999.

8. A. A. Samarskii, P. P. Matus, and P. N. Vabishchevich, Difference Schemes With Operator Factors, Kluwer, Dordreht, 2002.

9. A. A. Samarskii and P. N. Vabishchevich, Additive Schemes for Problems of Mathematical Physics, Moscow: Nauka, 1999 (in Russian).

10. R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, American Mathematical Society, 2001.

11. A. Toselli and O. B. Widlund, Domain Decomposition Methods—Algorithms and Theory, Springer-Verlag, 2005.

12. P. N. Vabishchevich, Additive schemes for certain operator-differential equations, Comput. Math. Math. Phys. 50(12) (2010), 2033–2043.

Keldysh Institute of Applied Mathematics (Received 03 02 2011) 4 Miusskaya Square 125047 Moscow Russia

 vabishchevich@gmail.com