- Email: [email protected]
Optimal error estimate for a mixed finite element method for compressible Navier–Stokes system
Applied Numerical Mathematics 45 (2003) 275–292 www.elsevier.com/locate/apnum
Optimal error estimate for a mixed finite element method for compressible Navier–Stokes system ✩ Jae Ryong Kweon Department of Mathematics, Pohang University of Science and Technology, Pohang 790-784, South Korea
Abstract A linearized stationary compressible viscous Navier–Stokes system is considered. A mixed finite element method is applied and the unique existence of the solution is established by the inf–sup condition. The convection terms, especially in the continuity equation, were thought of causing non-optimal order convergence, but in this paper error estimates of optimal order are derived by implementing the lowest order Raviart–Thomas elements. The error estimates for the normal and tangential components of velocity are also optimal on the interfaces of the interior triangles. It turns out that the non-symmetric discrete system can be reformulated into a symmetric form. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Compressible flows; Mixed finite element; Optimal convergence
1. Introduction So far we have not yet given a complete numerical understanding of the boundary value problems for the Navier–Stokes or the Stokes system for compressible flows [19,21]. Among many unsolved numerical problems is how to construct a finite element solution having optimal error estimate when a finite element method is applied to the system. Note that the system is an elliptic–hyperbolic type. For this, in this paper a mixed finite element method that considers the lowest order Raviart–Thomas elements [8] is applied to a linearized compressible viscous Navier–Stokes equations which can be obtained by linearizing the barotropic stationary compressible viscous Navier–Stokes equations around an ambient flow (U , P ). The equations are −µu − ν∇ div u + ρ (U · ∇)u + ∇p = f in Ω, (1.1) div u + κU · ∇p = g in Ω, u=0 on Γ, ✩
This work was supported by Korea Research Foundation Grant(KRF-2000-015-DP0012). E-mail address: [email protected] (J.R. Kweon).
0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00213-1
276
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
where Ω is a bounded domain with boundary Γ in R2 ; u = [u1 , u2 ] is a velocity vector, p is the pressure; ρ = ρ(P ) is a positive strictly increasing function that provides density as a function of pressure; ρ = dρ/dP , κ = ρ /ρ; the functions U = [U, V ], with U |Γ = 0, and P are the ambient flow which
the functions f and g are given functions, and the numbers µ and ν are the are C 1 functions in Ω; viscous constants with µ > 0 and µ + ν > 0. It is assumed that ρ and κ are given nonzero constants. Otherwise, system (1.1) becomes an incompressible Stokes flow. For the derivation of the equations in (1.1), we refer to [6,13,17]. Note that (1.1) may be obtained by simply adding additional convection terms in the momentum and continuity equations of the (incompressible) Stokes ones. It is often referred to as compressible Navier–Stokes equations. Also note that (1.1) is neither elliptic nor hyperbolic but contains features of each class of equations. In recent years, mixed discretizations have been used in many applications, in particular, for such problems where a major interest lies in the gradient variable rather than the primal one, for example, the flux in stationary incompressible flow problems or neutron diffusion or elasticity problems [2,4,20]. Recently in [16] a mixed scheme was applied to a simple compressible Stokes equation, so an error estimate was shown but not optimal. The convection terms in the equations were thought of causing not optimal order convergence. However, in this paper we show an optimal order error estimate for the solution of the mixed scheme of (1.1) by implementing the lowest order Raviart–Thomas finite elements. In [13] error estimates for the system were given, using an abstract formulation but not optimal. In [15] a weak formulation is proposed for the system with inflow boundary condition and an optimal error estimate is obtained. In [14] a discontinuous Galerkin method was applied to the system. For the stable finite element method, we refer to [3,12,5,9]. It has been observed in [10] that the finite volume methods are very similar to the Raviart–Thomas finite element when they were applied to the Poisson equation. Later in [11] this observation was generalized to the Stokes and elasticity problems by introducing a (new) finite element scheme. This scheme exhibits local properties of conservation of mass and momentum: In fact the method may be viewed as a finite volume scheme. In this paper, based on the integration by part formula (1.3) we formulate a mixed variational scheme and show stability results for the continuous and discrete problems, and also derive optimal error estimates in the variables of the velocity, the gradient of the velocity, and the pressure. Furthermore it is shown that the resulted non-symmetric discrete form for the mixed scheme of (1.1) can be cast into a symmetric matrix form (see (3.62) and (3.64)). Setting the gradient of the velocity u = (u1 , u2 ) by a tensor variable σ , system (1.1) is equivalently written as σ = ∇u in Ω, − div µ σ + ν tr(σ ) − p δ + ρ(U · ∇)u = f in Ω, (1.2) in Ω, tr(σ ) + κU · ∇p = g u=0 on Γ, where σ = (σij )1i,j 2 , σij = ∂ui /∂xj , δ is the identity tensor, and for a tensor τ = (τij )1i,j 2 ,
(div τ )i =
2 ∂τij j =1
∂xj
,
i = 1, 2,
and tr(τ ) = τ11 + τ22 is the trace of the tensor τ .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
277
In formulating a weak formulation of (1.2), the following integration by part formula will be used: σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx = 0. (1.3) Ω
Ω
For showing the stability result, some restrictions are needed. Let µ0 = min{µ, µ + 2ν}, µ1 = µ + 2|ν| and γ∗ = | div U |∞ max{ρ, κ}. In this paper we assume that the numbers µ0 , µ1 and γ∗ satisfy cγ∗ <
µ0 µ21
and
1 C < , µ0 γ∗
(1.4)
where c and C are generic constants. Note that when ν = 0, (1.4) becomes cγ∗ < 1/µ < C/γ∗ . In this paper a finite element scheme for (1.2) is given and unique existences of the continuous and discrete solutions are shown. Our main result is, which is derived in Section 3. Let (σ˜ = (σ, p), u˜ = (u, λn , λt )) and (σ˜ h = (σh , ph ), u˜ h = (uh , λnh , λt h )) be the solutions of (2.11) and (3.1), respectively. Define
κ∗ = Cµ0 µ−2 α∗ = µ(µ + 2ν)/ 3(µ + ν)2 + 5ν 2 . µ∗ = µ0 − Cγ∗ , 1 − γ∗ , If (1.4) holds and µ + 2ν = 0, then the following error estimates hold: ˆ 0,Ω , µ∗ σ − σh 0,Ω + κ∗ p − ph 0,Ω C σ − σˆ 0,Ω + p − p ˆ 0,Ω ), α∗ u − uh 0,Ω C(σ − σh 0,Ω + u − u
(1.5)
where (σˆ , p) ˆ ∈ Xh and uˆ ∈ V 1h are the interpolations of (σ, p) and u, defined by (3.9) and (3.12), respectively. Note that the numbers µ∗ and κ∗ are well-defined under (1.4), and also α∗ = 0 with µ + 2ν = 0. Note that (1.5) is optimal in a sense that the errors in both sides were measured with the same norms. We give the error estimate for the normal and tangential components of velocity on the interelement boundaries. Let λn = u · n and λt = u · t be the normal and tangential components of u on the boundary ∂K of the triangle K ∈ Th where Th is a given triangulation. Assume that µ + 2ν = 0 and (1.4) holds. If div U is constant on each triangle K ∈ Th , then 0 ρ (u · t) − λt h + ρh0 (u · n) − λnh 1/q,p,∂K Ch2/p Z1 , (1.6) h 1/q,p,∂K where Z1 = |σ |1,Ω + |p|1,Ω + | div(µ σ + ν tr(σ )δ)|1,Ω and ρh0 is the projection operator defined in (3.14). For k 0 and 1 p < ∞ we denote by W k,p (Ω) the Sobolev space with norm wk,p,Ω and seminorm |w|k,p,Ω (see [1]). Let W 1/q,p (Γ ) be the space of the traces v|Γ over Γ of the functions 1,p v ∈ W 1,p (Ω) where 1/p + 1/q = 1. Also W0 (Ω) denote the space of functions v in W 1,p with 1,p zero boundary value. W −1,q (Ω) and W −1/q,q (Γ ) denote the dual spaces of W0 (Ω) and W 1/q,p (Γ ), 1,2 respectively. If p = 2, W k,2 (Ω) and W0 (Ω) are denoted by H k (Ω) and H01 (Ω), respectively, and the norm uk,2,Ω is replaced by uk,Ω . Also L20 (Ω) denote the space of functions in L2 (Ω) with zero mean value over Ω. We denote by
2 Wq (div; Ω) = v ∈ Lq (Ω) ; div v ∈ Lq (Ω) , where 1 q < ∞, and W2 (div; Ω) means H(div; Ω). Finally we set
Mβ = χ ∈ L2 (Ω): χβ,Ω < ∞ ,
278
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
with norm χβ,Ω = χ20,Ω + div(U χ)20,Ω . We use the Sobolev embedding theorems: W 1,p (D) 1→ Lq (D) (p q 2p/(2 − p)), and W 1/q,p (D) 1→ L2 (∂D) (4/3 < p < q), where D ⊂ R 2 is a bounded domain. For the tensor τ = (τij )1i,j 2 we denote by τi , i = 1, 2, the vector τi = (τi1 , τi2 ) and the normal trace τ n = (τ1 · n, τ2 · n). The normal and tangential components, τnn and τnt , of the vector τ n shall be defined by τnn = τ n · n =
2
τ1j nj n1 +
j =1
τnt = τ n · t =
2 j =1
τ1j nj t1 +
2
τ2j nj n2 ,
j =1 2
τ2j nj t2 ,
j =1
where t denotes the unit tangential vector to Γ. The paper is organized as follows. In Section 2, based on the integration by part formula, several bilinear forms are defined and a variational formulation for (1.2) is presented, which consider the Lagrange multipliers on inter-element boundaries and an inf–sup condition is shown for the stability of the continuous problem. In Section 3 the continuous problem is discretized and the stability result is shown by the discrete inf–sup condition. Finally optimal error estimates are derived and the nonsymmetric matrix generated by the discrete problem is reformulated into a symmetric form. We use C to denote a generic positive constant, not depending on the mesh size h. Note that C may take different values in different places.
2. Mixed formulation In this section we first give a weak form for (1.2) (see (2.2)), where in [16] a similar form was analyzed, and consider the weak form that employs the Lagrange multipliers on the interelement boundaries (see (2.6)) and show a unique existence. We next give a regularity result for (1.1), which can be derived from [6]: Theorem A. Let Ω be a bounded domain of class C 2 and q > 2. Assume that f ∈ (Lq (Ω))2 and g ∈ W 1,q (Ω) with Ω g dx = 0. Then there is a constant K1 > 0 such that if U 2,q,Ω < K2 for a constant K2 , there is a unique solution (u, p) ∈ (W 1,q (Ω))2 × W 0,q (Ω) of (1.1), satisfying u2,q,Ω + p1,q,Ω K1 (f 0,q,Ω + g1,q,Ω ). For the weak formulation of (1.2), we define Qβ = Mβ ∩ L20 (Ω), V β = (Mβ )2 and
2 X0 = (τ, χ) ∈ (V 0 )2 × Q0 : µτ + ν tr(τ ) − χ δ ∈ H (div; Ω) .
(2.1)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
279
Note that V0 = L2 (Ω) and Q0 = L20 (Ω). Using (1.3) and (2.1), a weak formulation of (1.2) reads: find ((σ, p), u) ∈ X0 × V 0 such that σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx = 0, Ω Ω (2.2) div µσ + ν tr(σ ) − p δ · v + ρu · div(U v) dx = − f · v dx, Ω Ω tr(σ )χ − κp div(U χ) dx = gχ dx, Ω
Ω
for all (τ, χ) ∈ X0 , v ∈ V β and χ ∈ Qβ . It is obvious that (2.2) is equivalent to (1.2). In this paper we consider a finite element discretization for (2.2) by introducing the Lagrange multipliers on the interelement boundaries so that the normal trace (µτ + (ν tr(τ ) − χ)δ)n is continuous
This Lagrange multiplier is nothing else than the across interelement edges of the triangulation Th of Ω. normal component u · n and the tangential component u · t of the velocity. The tensor variable will be approximated by the lowest degree of the Raviart–Thomas finite element [8] 2 RT0 (K) := P0 (K) + x · P0 (K), where P0 (K) is the space consisting of constant functions on K and x = (x1 , x2 ) ∈ R2 , while the velocity and the pressure will be constant on each element of the triangulation and the Lagrange multiplier will be constant on every edge of each element. Note that the space (H (div; K))2 , K ∈ Th , is not appropriate in obtaining an approximation of the tensor τ by the Raviart–Thomas finite elements. In fact, if an interpolation operator is to be defined on the space (RT0 (K))2 by using the three normal components of a tensor τ on ∂K, K ∈ Th as the degrees of freedom, then the tensor τ should be slightly smoother than merely belonging to (H (div; K))2 [8]. From this point of view we will require q > 2 in the index q of the space (Wq (div; K))2 , which can be considered a finite element space for the tensor variable. Now we will choose the corresponding functional spaces that the Lagrange multipliers can be defined on interelement edges of the triangulation Th and the normal trace (µτ + (ν tr(τ ) − χ)δ)n is continuous there. Let p be a real number with 4/3 < p < 2 and q its conjugate where 1/p + 1/q = 1. We consider the following spaces [11]: Let V 1 = (Lq (Ω))2 , M = L20 (Ω) ∩ Lq (Ω) and
2 X = (τ, χ) ∈ V 1 × M: µτ + ν tr(τ ) − χ δ |K ∈Wq (div; K), ∀K ∈ Th , (2.3) = η ∈ W 1/q,p (∂K): η|∂K∩Γ = 0, ∀K ∈ Th . V 2 K∈Th
1 = V β ∩ (Lq (Ω))2 and We also denote by V V = V 1 × (V2 )2 ,
=V 1 × (V2 )2 , V
Q = Qβ ∩ Lq (Ω),
(2.4)
where Vβ and Qβ are defined above. The following variables are introduced: σ˜ = (σ, p),
τ˜ = (τ, χ),
u˜ = (u, λn , λt ),
v˜ = (v, ηn , ηt ).
(2.5)
280
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using the spaces in (2.3), (2.4) and the variables given in (2.5), a weak formulation for (1.2) reads: find σ˜ = (σ, p) ∈ X, u˜ = (u, λn , λt ) ∈ V such that σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx Ω K∈Th K τnt λt ds = 0, ∀τ˜ ∈ X, µτ + ν tr(τ ) − χ δ nn λn ds − µ − ∂K ∂K K∈ T K∈ T h h 1, div µσ + ν tr(σ ) − p δ · v + ρu div(U v) dx = − f v dx, ∀v ∈ V K Ω K∈Th (2.6) η ds = 0, ∀η ∈ V , µσ + ν tr(σ ) − p δ n n 2 nn K∈Th ∂K σnt ηt ds = 0, ∀ηt ∈ V2 , µ ∂K K∈ T h χ tr(σ ) − κp div(U χ) dx = gχ dx, ∀χ ∈ Q. K Ω K∈Th
Note that the Lagrange multipliers were introduced in (2.6) but not in (2.2). In order to write (2.6) in a compact form we consider four bilinear forms on the functional spaces in (2.3) and (2.4): (2.7) a(σ˜ , τ˜ ) = σ : µτ + ν tr(τ ) − χ δ dx, σ˜ , τ˜ ∈ X, Ω
and b(τ˜ , v˜ ) =
K∈Th K
−µ
div µτ + ν tr(τ ) − χ δ · v dx − µτ + ν tr(τ ) − χ δ nn ηn ds
K∈Th∂K
τnt ηt ds,
τ˜ ∈ X, v˜ ∈ V ,
(2.8)
K∈Th∂K
and c(p, χ) = −
p div(U χ) dx,
p ∈ M, χ ∈ Q,
(2.9)
K∈Th K
and d(σ, χ) =
χ tr(σ ) dx,
σ ∈ (V 1 )2 , χ ∈ M.
(2.10)
K∈Th K
Using (2.7)–(2.10), (2.6) can be written by: find σ˜ ∈ X and u˜ ∈ V such that ˜ = 0, τ˜ ∈ X, a(σ˜ , τ˜ ) + b(τ˜ , u) , b(σ˜ , v˜ ) − c(ρ u, v) = −f , v , ∀v˜ ∈ V c(κp, χ) + d(σ, χ) = g, χ , ∀χ ∈ Q, where u, v = Ω uv dx.
(2.11)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
281
Note that the form a(σ˜ , τ˜ ) is not symmetric on X × X. However, this non-symmetric form a can be decomposed into a symmetric form plus the other bilinear one by taking the test function τ˜ ∈ X suitably. For details, see (3.58). In order to show the unique existence of solution of (2.11), which is shown in Theorem 2.1 below, we recall the inf–sup condition for the bilinear form b, which can be shown by the same argument as given in [16, Lemma 2.1]: There is a positive constant α such that sup τ˜ ∈X τ˜ =0
b(τ˜ , v˜ ) α v0,Ω , τ˜ X
∀v˜ ∈ V ,
(2.12)
where τ˜ X := (τ 20,Ω + χ20,Ω )1/2 . Now using (2.12) and Theorem A, we show a unique existence of the solution for the continuous problem (2.11): Theorem 2.1. Assume that the conditions given in Theorem A hold. Assume that β2,q,Ω is small in a suitable sense. Let (u, p) be the solution of (1.1) and σ = ∇u. Then (σ˜ = (σ, p), u˜ = (u, u · n, u · t)) is the unique solution of (2.11). Furthermore, if (1.4) holds, then µ∗ σ 0,Ω + κpβ 0,Ω + κ∗ p0,Ω C f 0,Ω + g0,Ω , √ (2.13) αu0,Ω 2(µ1 + 1)σ 0,Ω , √ where µ∗ = µ0 − Cαγ∗ and κ∗ = Cµ0 µ−2 2 − γ∗ with µ2 = µ1 + ρ|U |∞ . Proof. Since the conditions given in Theorem A is assumed, the pair σ˜ = (σ, p) ∈ X, u˜ = (u, u · n, u · t) ∈ V and is the solution of (2.11). The uniqueness follows from Theorem A. The inequality (2.13), constants µ∗ and κ∗ can be derived by applying the same procedures as given in the proof of [16, Theorem 2.1]. ✷
3. Mixed discretization Before proceeding further, we recall that in [16] a mixed finite element scheme that implements the Raviart–Thomas elements of order l was applied to a (simple) compressible Stokes system but does not consider the Lagrange multipliers on the interelement boundaries like (2.6). However, since the governing equations and finite elements considered in this paper have generically similar characters like [16], we often use some similar arguments as employed in [16] due to technicality and complexity. Now we are going to approximate the solution of (2.11) by implementing the lowest order Raviart– Thomas finite elements. Let h be an approximate indicator. With finite-dimensional spaces Xh ⊂ X, V h ⊂ V , Mh ⊂ M and Qh ⊂ Q, we find the approximate solutions σ˜ h ∈ X h and u˜ h ∈ V h such that ∀τ˜h ∈ X h , a(σ˜ h , τ˜h ) + b(τ˜h , u˜ h ) = 0, (3.1) b(σ˜ h , v˜ h ) − c(ρ uh , vh ) = −f , vh , ∀v˜ h ∈ V h , ∀χh ∈ Qh . c(κph , χh ) + d(σ˜ h , χh ) = g, χh , In order to analyze the discrete problem (3.1) we are going to specify the finite-dimensional spaces Xh , V h , Mh , Qh and their degrees of freedom. Let Th be a given regular triangulation for the closure of Ω and Pl (K) the space of polynomials of degree l, where K ∈ Th is a triangle.
282
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Let RT0 (K) denote by the lowest order Raviart–Thomas element 2 RT0 (K) := P0 (K) + x · P0 (K). Note that any v ∈ RT0 (K) is uniquely determined by its normal components n · v |ei on the edges ei , 1 i 3, of K ∈ Th , where n denotes the outer normal vectors of K. Letting D := (RT0 (K))2 , we define an operator ΠK1 : (Wq (div; K))2 → D by ΠK1 τ ∈ (RT0 (K))2 such that 2 1 (3.2) ΠK τ − τ n · v 0 ds = 0, ∀v 0 ∈ R0 (∂K) , ∂K
where R0 (∂K) = {v: v|e ∈ P0 (e), ∀e ∈ ∂K}. Note that the boundary integral in (3.2) may be written by 1 (3.3) ΠK τ n · v 0 ds = τ n, v0 ∂K , ∂K
where . , . ∂K denotes the duality pairing between (W −1/q,q (∂K))2 and (W 1/q,p (∂K))2 . The interpolation error estimates for ΠK1 are cited in [11]: if Th is regular, then 1 Π τ − τ ChK |τ |1,K , K 0,K 1 div Π τ − τ C div τ 0,K , K
0,K
(3.4) (3.5)
where hK is a diameter of K. Using the operator ΠK1 of (3.3), let Πh1 be defined by Πh1 τ|K = ΠK1 τ,
∀K ∈ Th ,
(3.6)
for all τ ∈ {τ ∈ (L (Ω)) : τ ∈ (Wq (div; K)) , ∀K ∈ Th }. Define an approximate space of X as follows:
(3.7) X h = τ˜h ∈ X: τh|K ∈ D, χh|K ∈ P0 (K), ∀K ∈ Th , q
4
2
where τ˜h = (τh , χh ). Let ρh be an projection operator defined by 1 rh χ dx, ρh χ = rh χ − |Ω|
(3.8)
Ω
where rh is the L2 projection onto the space consisting of piecewise constants. Noting that (µ τh + (ν tr(τh ) − χh )δ)|K ∈ (RT0 (K))2 for ∀K ∈ Th , we define an operator Πh : X → X h , Πh τ˜ = (τh , χh ) ∈ X h (see [11]) such that χh = ρh χ and for ∀v ∈ (R0 (∂K))2 and ∀K ∈ Th ,
µτ + ν tr(τ ) − χ δ − µτh + ν tr(τh ) − χh δ n · v ds = 0. (3.9) ∂K
Using (3.7) and (3.9) we have µτh + ν tr(τh ) − χh δ = Πh1 µτ + ν tr(τ ) − χ δ , and the interpolation error estimate for Πh is given in [11]: assuming that Th is regular and (τ, χ) ∈ (H 1 (Ω))4 × H 1 (Ω), we have (3.10) Πh τ˜ − τ˜ 0,Ω Ch |τ |1,Ω + |χ|1,Ω .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
For the approximations of V , M and Q, we set
2 V 1h = v h ∈ V 1 : v h |K ∈ P0 (K) , ∀K ∈ Th ,
V2h = ηh ∈ V2 : ηh |∂K ∈ R0 (∂K), ∀K ∈ Th , V h = V 1h × (V2h )2 , M = χ ∈ M: χ | ∈ P (K), ∀K ∈ T }, Qh = Mh ∩ Q. h h K 0 h Define the operator Ph1 : V 1 → V 1h by 2 Ph1 v|K = PK1 v, ∀v ∈ V 1 , PK1 v ∈ P0 (K) such that 2 (PK1 v − v) · v 0 dx = 0, ∀v 0 ∈ P0 (K) , ∀K ∈ Th .
283
(3.11)
(3.12)
K
Consider the operator Ph : V → V h defined by Ph v˜ = Ph (v, ηn , ηt ) = Ph1 v, ρh0 ηn , ρh0 ηt , where ρh0 is the projection operator from K∈Th L2 (∂K) onto K∈Th R0 (∂K), i.e., 0 0 ρh : V2 → V2h , ρh η ∈ R0 (∂K), 0 ρh η − η v0 ds = 0, ∀v0 ∈ R0 (∂K), ∀K ∈ Th .
(3.13)
(3.14)
∂K
Note that from the Sobolev’s imbedding theorem W 1/q,p (∂K) 1→ L2 (∂K)(4/3 < p < q), the operator ρh0 is well-defined on V2 (cf. see [11]). In next lemma we establish a discrete inf–sup condition for the bilinear form b, in which the proof is similar to the ones in [11] and [16, Lemma 3.1]: Lemma 3.1. Assume that Th is regular and µ + 2ν = 0. Then there exists a constant α∗ independent of h such that b(τ˜h , v˜ h ) α∗ v h 0,Ω , ∀v˜ h ∈ V h , (3.15) sup τ˜h ∈X h τ˜h X −1/2
where α∗ = C/µ3
with µ3 defined in (3.16).
Proof. Let v˜ h = (v h , µnh , µt h ) ∈ V h . By [11, Proposition 3.1] there is a τh∗ ∈ (Wq (div; Ω))2 such that τh∗ |K ∈ (RT0 (K))2 , ∀K ∈ Th , div τh∗ = v h in Ω and τh∗ 0,Ω Cv h 0,Ω . Picking a function τh satisfying µτh + ν tr(τh )δ = τh∗ , the components of τh can be expressed in terms of those of τh∗ and (µ + ν)2 + ν 2 ∗ 2 ∗ 2 + τ τ h11 0,Ω h22 0,Ω (µ2 + 2µν)2 1 ∗ 2 ∗ 0,Ω + τh21 20,Ω µ3 τh∗ 20,Ω Cµ3 v h 20,Ω , + 2 τh12 µ
τh 20,Ω 2
(3.16)
where µ3 = (3(µ + ν)2 + 5ν 2 )/(µ2 (µ + 2ν)2 ). Finally, letting τ˜h = (τh , 0) ∈ X h , we have b(τ˜h , v˜ h ) = −1/2 v h 20,Ω (since τh ∈ (Wq (div; Ω))2 ). Thus the required inequality (3.15) follows with α∗ = Cµ3 . ✷
284
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using Lemma 3.1 and considering some similar arguments used in the proof of [16, Theorem 3.1], we show a unique existence of the discrete solution of (3.1). Theorem 3.1. Let µ + 2ν = 0. Assume that (1.4) holds. Then problem (3.1) has a unique solution (σ˜ h , u˜ h ) ∈ X h × V h , satisfying µ∗ σh 0,Ω + κ∗ ph 0,Ω C f 0,Ω + g0,Ω , (3.17) α∗ uh 0,Ω Cσh 0,Ω , where µ∗ and κ∗ are given in (3.26), and α∗ is given in (3.15), and C is a constant not depending on h. Proof. Let Zh = {τ˜h ∈ X h : b(τ˜h , v˜ h ) − c(ρuh , v h ) = −f , vh , ∀v˜ h ∈ V h } and let τ˜h = (τh , χh ) ∈ Zh be arbitrary. Recalling that V 1h consists of piecewise constant functions, and letting (div U )h and f h be the L2 projections into the piecewise constant functions which are constants on each triangle K ∈ Th , we have (3.18) div µτh + ν tr(τh ) − χh δ + ρuh (div U )h + f h v h dx = 0, K∈Th K
for all v˜ h = (vh , 0, 0) ∈ V h . It follows from (3.18) that div µτh + ν tr(τh ) − χh δ + ρuh (div U )h = −f h .
(3.19)
Setting σ1 = µτh + (ν tr(τh ) − χh )δ with χh ∈ L20 (Ω), we have div σ1 = −ρuh (div U )h − f h . We now need to show that (3.20) χh 0,Ω c1 µ1 τh 0,Ω + c2 γ∗ uh 0,Ω + f h 0,Ω , ∀τ˜h ∈ Zh , √ √ and c2 = C/ 2. Indeed, to show (3.20), we define where c1 = 2(2 + C) µ + 2ν |Ω|−1 tr(τh ) dx, L= 2 Ω
and σ2 = µτh + (ν tr(τh ) − χh − L)δ. Using (3.19), we have div σ2 = −ρuh (div U )h − f h , and Ω tr(σ2 ) dx = 0. But using [8, p. 161, Proposition 3.1], there is a constant C such that σ2D + div σ2 0,Ω , (3.21) σ2 0,Ω C 0,Ω where σ2D = σ2 − 12 tr(σ2 )δ. On the other hand, µ σ2D = µτh − tr(τh )δ, 2 and using the similar procedures as given in [16, Theorem 3.1], h 0,Ω + C div σ2 0,Ω + Lδ0,Ω σ1 0,Ω µCτ + 1 + 2|ν| τh 0,Ω + C γ∗ uh 0,Ω + f h 0,Ω . µ 2C Since χh δ = µτh + ν tr(τh )δ − σ1 , and using (3.23) we have 1 χh 0,Ω √ µ1 τh 0,Ω + σ1 0,Ω c1 µ1 τh 0,Ω + c2 γ∗ uh 0,Ω + f h 0,Ω , 2
(3.22)
(3.23)
(3.24)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
285
√ √ and c2 = C/ 2. Thus (3.20) follows. Now, recalling µ0 = min{µ, µ + 2ν} and where c1 = 2(2 + C) using (3.20), we obtain f , uh + g, ph = a(σ˜ h , σ˜ h ) + d(σh , ph ) + c(ρuh , uh ) + c(κph , ph ) µ0 σh 20,Ω − γ∗ uh 20,Ω + ph 20,Ω µ0 µ0 −1 2 − α∗ γ∗ σh 0,Ω + − γ ph 20,Ω 2 4(c1 µ1 )2 µ0 2 2 2 − u + f γ h h ∗ 0,Ω 0,Ω 2(c1 µ1 )2 1 (3.25) µ∗ σh 20,Ω + κ∗ ph 20,Ω − Cf h 20,Ω 2 where µ∗ = (µ0 − c3 γ∗ ) with c3 = 2α∗−1 + µ0 γ∗ /(c1 µ1 )2 , and κ∗ = µ0 /4(c1 µ1 )2 − γ∗ . Finally, computing f , uh + g, ph and using f h 0,Ω f 0,Ω , the inequality (3.17a) follows from (3.25). Also (3.17b) follows from Lemma 3.1 and Eq. (3.1a). The uniqueness follows from (3.17). ✷ From the definition of the operators Πh and Ph by (3.9) and (3.13), respectively, we recall that the bilinear form b defined by (2.8) satisfies b(Πh τ˜ − τ˜ , v˜ h ) = 0, b(τ˜h , Ph v˜ − v˜ ) = 0,
∀v˜ h ∈ V h , ∀τ˜ ∈ X, ∀v˜ ∈ V , ∀τ˜h ∈ X h .
(3.26)
We are now ready to derive an optimal error estimate which can be compared with [16, Theorems 3.2 and 4.2]. The derivation is similar to the one in [16, Theorems 3.2]. Theorem 3.2. Assume that (1.4) holds and µ + 2ν = 0. Let Th is a given regular triangulation of the domain Ω. Let (σ˜ , u) ˜ be the solution of (2.11) and (σ˜ h , u˜ h ) the solution of (3.1). Assume that 1 4 σ˜ = (σ, p) ∈ (H (Ω)) × (H 1 (Ω) ∩ L20 (Ω)). Then there exists a constant C independent of h such that ˆ 0,Ω + u − u ˆ 0,Ω , µ∗ σ − σh 0,Ω + κ∗ p − ph 0,Ω C σ − σˆ 0,Ω + p − p (3.27) ˆ 0,Ω , α∗ u − uh 0,Ω C σ − σh 0,Ω + u − u where σˆ , pˆ and uˆ are the interpolations of σ, p and u, respectively. Proof. From (2.11) and (3.1), and using (3.26), we have ˜ = a(σ˜ − Πh σ˜ , τ˜h ), ∀τ˜h ∈ X h , a(σ˜ h − Πh σ˜ , τ˜h ) + b(τ˜h , u˜ h − Ph u) 1 b(σ˜ h − Πh σ˜ , v˜ h ) − c uh − Ph u, vh = −c u − Ph1 u, v h , ∀v˜ h ∈ V h , c(ph − ρh p, χh ) + d σh − Πh1 σ, χh = c(p − ρh p, χh ) + d σ − Πh1 σ, χh , Applying to (3.29) the same procedures used in showing (3.20), we have ˆ 0,Ω c1 µ1 σh − Π 1 σ 0,Ω + c2 γ∗ uh − P 1 u + u − P 1 u ph − p h
h
0,Ω
h
(3.28) (3.29) ∀χh ∈ Qh .
, 0,Ω
(3.30)
(3.31)
where ci (i = 1, 2) are given in (3.20). Taking v˜ h = u˜ h − Ph u˜ in (3.28), τ˜h = σ˜ h − Πh σ˜ in (3.29), and χh = ph − pˆ with pˆ = ρh p in (3.30), and adding the resulted equations, and applying the same procedures used in showing (3.25), and using (3.31), we obtain
286
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
a(σ˜ h − Πh σ˜ , σ˜ h − Πh σ˜ ) + d σh − Πh1 σ, ph − pˆ ˆ ph − p) ˆ + c uh − Ph1 u, uh − Ph1 u + c(ph − p, 2 2 1 ˆ 20,Ω − C u − Ph1 u 0,Ω , (3.32) µ∗ σh − Πh1 σ 0,Ω + κ∗ ph − p 2 where µ∗ and κ∗ are defined in (3.25). On the other hand, the resulting right-hand is estimated by a(σ˜ − Πh σ˜ , σ˜ h − Πh σ˜ ) + c u − Ph1 u, uh − Ph1 u ˆ + d σ − Πh1 σ, ph − pˆ + c(p − p, ˆ ph − p) 2 2 ˆ 20,Ω + ε3 uh − Ph1 u 0,Ω ε1 σh − Πh1 σ 0,Ω + ε2 ph − p 2 2 + C σ − Πh1 σ 0,Ω + p − p ˆ 20,Ω + u − Ph1 u 0,Ω , ∀εi > 0, i = 1, 3, (3.33) where C = C(γ∗ , µ1 , |U|∞ ). Furthermore, using Lemma 3.1 and (3.28), , α∗ uh − P 1 u C σh − Π 1 σ + u − P 1 u h
h
0,Ω
0,Ω
h
0,Ω
(3.34)
where α∗ is given in (3.15). Thus, combining (3.32)–(3.34) and using the triangle inequality, the required inequality follows. ✷ In the formulation of (2.11), λn (respectively, λt ) is nothing else than the normal (respectively, the tangential) component u · n (respectively u · t) of the velocity. Therefore, in the discrete problem (3.1), λnh (respectively, λt h ) is an approximation of u · n (respectively, u · t). We are now going to derive error bounds for λnh − u · n and λt h − u · t. To do this, we cite the following lemma (see [11]): Lemma 3.2. For any T ∈ (W −1/q,q (∂K))2 , there exists a unique τh ∈ (RT0 (K))2 such that 2 (T − τh n) · v0 ds = 0, ∀v0 ∈ R0 (∂K) , ∂K
where R0 (∂K) = {v: v|e ∈ P0 (e), ∀e ∈ ∂K}. Furthermore, assume that Th is uniformly regular. Then there is a constant C independent of h such that τh 0,q,K CT −1/q,q,∂K , τh 0,K Ch2/p−1T −1/q,q,∂K , div τh 0,K Ch2/p−2T −1/q,q,∂K , where p and q are conjugate. Using Lemma 3.2 and applying exactly the same arguments given in [11, Theorem 3.3], we obtain: Theorem 3.3. Let (σ, u) be the solution of (2.11) and (σ˜ h , u˜ h ) the solution of (3.1). Assume that µ + 2ν = 0. If Th is uniformly regular and σ˜ = (σ, p) ∈ (H 1 (Ω))4 × (H 1 (Ω) ∩ L20 (Ω)), then there exists a constant C independent of h such that 0 ρ (u · t) − λt h + ρ 0 (u · n) − λnh Ch2/p−1 Z, (3.35) h
1/q,p,∂K
h
for all K ∈ Th , where Z = |σ |1,Ω + |p|1,Ω .
1/q,p,∂K
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
287
The error estimate (3.35) is not optimal. In order to get an optimal estimate we need the following Lemmas: Lemma 3.3. (1) Assume that |∇ U |∞ is small enough. Then there is a unique solution z of the problem: µz + ν∇ div z + ρ div(U z) = w z=0
in Ω, on Γ,
(3.36)
satisfying µz1,Ω Cw0,Ω , with C = C(|∇U |∞ ). Also the solution z satisfies µzk,q,Ω Cwk−2,q,Ω with C = C(Ω, U , k), k is an integer 1, 1 < q < ∞. (2) Let k −1 be an integer and let q ∈ (2/(k + 2), ∞). Assume that U ∈ W k+3,q (Ω) with U |Γ = 0 and Γ ∈ C k+3 . Assume that µ + 2ν = 0 and 2/(µ + 2ν) > αk ≡ CU k+3,q,Ω . For any S ∈ W k,q (Ω) there is a solution s ∈ W k,q (Ω) of problem 2 s − κ div(U s) = S, µ + 2ν
(3.37)
2 − αk )sk,q,Ω CSk,q,Ω for a constant C. satisfying ( µ+2ν
Proof. (1) It easily follows from the standard Lq theory of elliptic equations and [17, Lemma 2.4]. (2) It follows from [7, Theorem 1.1] that problem (3.37) has a solution s satisfying the required inequality. For further references one may refer to Beirão da Veiga [6, p. 234] and Lax & Phillips [18]. ✷ On the basis of (3.36) and (3.37) we consider the following problem µz + ν∇ div z + ρ div(U z) = w in Ω, 2 s − κ div(U s) + div z = h in Ω, µ+2ν z=0 on Γ, s dx = 0,
(3.38)
Ω
where z and s are the unknown variables, and w, h are given functions. From Lemma 3.3 one can derive that if w ∈ L2 (Ω) and h ∈ H 1 (Ω) are given, then there is a unique solution (z, s) ∈ (H01 (Ω))2 × L20 (Ω) of (3.38) satisfying z2,Ω K1 w0,Ω , s1,Ω K2 h1,Ω + w0,Ω , where K1 is a constant and K2 = C(2/(µ + 2ν) − α1 )−1 . In order to derive a dual formulation of (1.2), let γ = (γij )1i,j 2 be a tensor variable and let µγ + ν tr(γ ) − s δ = µ∇z + ν div zδ,
(3.39) (3.40)
(3.41)
where δ is the identity tensor. Taking the trace to both sides of (3.41), we have tr(γ ) =
2 s + div z. µ + 2ν
(3.42)
288
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using (3.41) and (3.42), system (3.38) becomes µγ + ν tr(γ ) − s δ = µ∇z + ν div zδ div µγ + ν tr(γ ) − s δ + ρ div(U z) = w tr(γ ) − κ div(U s) = h z=0 s dx = 0.
in Ω, in Ω, in Ω, on Γ,
(3.43)
Ω
On the other hand, multiplying τ to both sides of (3.41) and doing the integration by part and using z|Γ = 0, we have the following formula (see (1.3)): τ : µγ + ν tr(γ ) − s δ dx + div µτ + ν tr(τ )δ · z dx = 0. Ω
Ω
Hence a dual formulation for (2.11) is: Find γ˜ = (γ , s) ∈ X, z˜ = (z, αn , αt ) ∈ V such that τ˜ ∈ X, a(τ˜ , γ˜ ) + b1 (τ, z˜ ) = 0, b(γ˜ , v˜ ) − c(ρv, z) = (w, v), ∀v˜ ∈ V , c(κ χ, s) + d(γ , χ) = (h, χ), ∀χ ∈ M, where div µτ + ν tr(τ )δ · z dx − µτ + ν tr(τ )δ nn αn + µτnt αt ds b1 (τ, z˜ ) = K∈Th K
(3.44)
K∈Th∂K
for τ˜ = (τ, χ) ∈ X and z˜ ∈ V ,
(3.45)
and the bilinear forms a, b, c and d are defined in (2.7)–(2.10), respectively. Using (3.39), (3.40) and (3.43c) we obtain the following lemma: Lemma 3.4. Let µ + 2ν = 0. Assume that w ∈ (L2 (Ω))2 and h ∈ H 1 (Ω). Assume that the solution (z, s) of (3.38) has the regularity µγ + (ν tr(γ ) − s)δ ∈ (H (div; Ω))2 and s ∈ L20 (Ω), where µγ + (ν tr(γ ) − s)δ = µ∇z + ν div zδ. Then γ˜ = (γ , s) and z˜ = (z, z · n, z · t) is the unique solution of (3.44) and also satisfy z2,Ω Cw0,Ω , γ 1,Ω + s1,Ω C h1,Ω + w0,Ω ,
(3.46) (3.47)
where C = C(K1 , K2 ). Using the estimates (3.46) and (3.47), optimal error estimates for u · n and u · t are derived: Theorem 3.4. Assume that the conditions of Theorem 3.3 hold and div[µσ + ν tr(σ )δ] ∈ (H 1 (Ω))2 . Assume that µ + 2ν = 0 and (1.4) hold. Then if div U is constant on each triangle K ∈ Th , there exists a constant C not depending on h such that 0 ρ (u · t) − λt h + ρ 0 (u · n) − λnh Ch2/p Z1 , h
1/q,p,∂K
h
1/q,p,∂K
where Z1 = |σ |1,Ω + |p|1,Ω + | div[µσ + ν tr(σ )δ]|1,Ω .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Proof. Let γ˜ = (γ , s) ∈ X, z˜ = (z, αn , αt ) ∈ V be the solution of ∀τ˜ = (τ, χ) ∈ X, a(τ˜ , γ˜ ) + b1 (τ, z˜ ) = 0, b(γ˜ , v˜ ) − c(ρ v, z) = (u¯ h − uh , v), ∀v˜ ∈ V , c(κ χ, s) + (tr(γ ), χ) = 0, ∀χ ∈ M
289
(3.48)
where u¯ h = Ph1 u (Ph1 is defined by (3.12)). Using u¯ h − uh ∈ (L2 (Ω))2 and Lemma 3.4, we have z1,Ω Cu¯ h − uh 0,Ω
and
γ˜ 1,Ω Cu¯ h − uh 0,Ω ,
(3.49)
where C is a positive constant not depending on h. First, noting that c(u¯ h − uh , z) = 0 (by the integration by parts and U |Γ = 0), using (3.48b) and (3.26), we have u¯ h − uh 20,Ω = b(Πh γ˜ , u˜ − u˜ h )
(using (3.28))
= a(σ˜ h − σ˜ , Πh γ˜ ) = a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + a(σ˜ h − σ˜ , γ˜ )
(using (3.48a))
= a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + b1 (σ − σh , z˜ ).
(3.50)
On the other hand, letting Ph z˜ = (¯zh , αhn , αht ), where z¯ h := Ph1 z and αh := z¯ h|∂K , we have b(Πh σ˜ − σ˜ h , Ph z˜ ) = b1 Πh1 σ − σh , Ph z˜ = b1 (σ − σh , Ph z˜ ). Hence it follows from (3.29) and (3.26) that b1 (σ − σh , Ph z˜ ) − ρc u − uh , Ph1 z = 0.
(3.51)
If div U is constant on each triangle K ∈ Th , then c(u − uh , Ph1 z) = c(u¯ h − uh , Ph1 z). Subtracting (3.51) from (3.50), we have u¯ h − uh 20,Ω = a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + b1 (σ − σh , z˜ − Ph z˜ ) + ρ c(u¯ h − uh , Ph1 z) (using (3.49)) Ch2 Z1 |γ˜ |1,Ω + |z|1,Ω + C| div U |∞ u¯ h − uh 20,Ω ,
(3.52)
where Z1 = |σ |1,Ω + |p|1,Ω + | div(µ σ + ν tr(σ )δ)|1,Ω and for (3.52) we used div µσ + ν tr(σ )δ − Πh1 µσ + ν tr(σ )δ · z − Ph1 z dx, b1 (σ − σh , z˜ − Ph z˜ ) = K∈Th K
which follows from the properties of the operator Ph and the fact µσ + ν tr(σ )δ ∈ (Wq (div; Ω))2 . Using (3.52) and (3.49) we have u¯ h − uh 0,Ω C1 h2 Z1 ,
(3.53)
where C1 = C/(1 − C| div U |∞ ). Using exactly the same procedures as ones given in [11, Theorem 3.4] and the inequality (3.38) given in [11, Theorem 3.3], we obtain the required inequality. ✷
290
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
As mentioned in Section 2 the form a(σ˜ , τ˜ ) is not symmetric. However, by taking the test function τ˜h ∈ X suitably, system (2.11) can be transformed into a symmetric system. For τ˜ = (τ, 0) ∈ X let ˜ = b(τ˜ , u) ˜ and for τ˜ = (0, χ) ∈ X let a2 (σ˜ , τ˜ ) = a(σ˜ , τ˜ ), b2 (τ˜ , u) ˜ = b(τ˜ , u). ˜ a1 (σ˜ , τ˜ ) = a(σ˜ , τ˜ ), b1 (τ˜ , u) Then a1 (σ˜ , τ˜ ) is a symmetric form for τ˜ = (τ, 0) ∈ X and a2 (σ˜ , τ˜ ) = −d(σ, χ) for τ˜ = (0, χ) ∈ X. Thus system (2.11) becomes: find σ˜ ∈ X and u˜ ∈ V such that ˜ = 0, τ˜ = (τ, 0) ∈ X, a1 (σ˜ , τ˜ ) + b1 (τ˜ , u) , (3.54) b(σ˜ , v˜ ) − c(ρu, v) = −f , v , ∀v˜ ∈ V ˜ = g, χ , ∀χ ∈ Q, c(κ p, χ) + b2 (τ˜ , u) where u, v = Ω uv dx. Similar for the discrete system (3.1). Indeed, letting τh = 0 in (3.1a) and recalling that Mh is the space consisting of piecewise constants, we get tr(σh )χh dx − χh λnh ds = 0, ∀χh ∈ Mh . (3.55) K∈Th∂K
Ω
Combining (3.55) and (3.1c), the discrete problem (3.1) can be reformulated as follows: µσ : τ + ν tr(σ ) tr(τ ) + div µτ + ν tr(τ )δ · uh h h h h h h Ω K K∈ T h − (τh )nt λt h ds = 0, µτh + ν tr(τh )δ nn λnh ds − µ ∂K ∂K K∈T K∈ T h h κ ph χh div U + χh λnh ds = gχh , − K Ω K∈ Th K∈Th ∂K div µσ + ν tr(σ )δ · v + ρ u · v div U = − f vh, h h h h h K Ω Ω K∈Th + ν tr(σ )δ η ds + ph ηnh ds = 0, µσ − h h nn nh ∂K ∂K K∈Th K∈Th (σh )nt ηt h ds = 0, −µ ∂K
(3.56)
K∈Th
for all (τ˜h = (τh , χh ); v˜ h = (v h , ηnh , ηt h )) ∈ X h × Mh . The algebraic equations generated by (3.56) are of the form σ 0 A 0 BT DT E T 0 HT 0 p G 0 C 0 0 u = F , B 0 Z D H 0 0 0 q 0 E 0 0 0 0 r 0
(3.57)
where σ, p, u, q and r denote by the degrees of freedom for σh , ph , uh , λnh and λt h , respectively. Since A is (symmetric) positive definite and block diagonal, the unknown σ can be eliminated, i.e., (3.58) σ = − A−1 B T u + A−1 D T q + A−1 E T r . Putting σ into the third and fourth equations of (3.57), we have
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
F = Z − BA−1 B T u − BA−1 D T q − BA−1 E T r, −1
−1
−1 T
DA B u + DA D q + DA E r − Hp = 0. T
T
291
(3.59) (3.60)
Set A∗ = Z − BA−1 B T . If | div U |∞ is small enough, then A∗ will be invertible and the unknown u is eliminated and −1 T −1 −1 T −1 u = A−1 ∗ BA D q + A∗ BA E r + A∗ F.
(3.61)
Putting (3.61) into (3.60), and inserting (3.61) into (3.58) and using the equation Eσ = 0 we have
q + Yr − Hp = F
1 , X T
2 , (3.62) Y q + Mr = F T H q − Cp = −G, where
−1 T T X = DA−1 B T A−1 , ∗ DA B −1 T
= X + DA D , X T −1 E + DA−1 E T , Y = D A−1 B T A−1 ∗ BA −1 T T + EA−1 E T , M = EA−1 B T A−1 ∗ EA B
1 = −DA−1 B T A−1 F ∗ F,
(3.63)
2 = −EA−1 B T A−1 F ∗ F.
in (3.63) is invertible, then (3.62) is written by If X
−1 H
−1 Y p Λ HT X −C + HT X 1 T −1 T = , T −1
r Λ M−Y X Y H X Y 2
(3.64)
2 − Y T X
−1 F
1 − g and Λ2 = F
−1 F
1 . where Λ1 = HT X
Conclusion The linearized steady-state barotropic compressible viscous Navier–Stokes equations (1.1) are a system of PDEs of mixed type; the momentum equations form an elliptic subsystem in velocity variable and the continuity equation is a hyperbolic equation in pressure. Hence it is a non-self adjoint boundary value problem. However, from the numerical approximation of the mixed finite element method (3.1) by the lowest order Raviart–Thomas finite elements we have seen that the resulted discrete equations in (3.56) generate a symmetric form (3.57), which can be transformed into (3.64) and also the approximate solution of the discrete problem (3.1) converges to the true solution optimally.
References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] D.N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: Implementation, postprocessing and error estimates, RAIRO Modél. Math. Anal. Numér. 19 (1985) 7–32.
292
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
[3] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344. [4] D.N. Arnold, R.S. Falk, A new mixed formulation for elasticity, Numer. Math. 16 (1971) 322–333. [5] I. Babuška, W.C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (1978) 736–754. [6] H. Beirão da Veiga, An Lp -theory for the n-dimensional, stationary, compressible Navier–Stokes equations, and the incompressible limit for compressible fluids. The equilibrium solutions, Commun. Math. Phys. 109 (1987) 229–248. [7] H. Beirão Da Veiga, Existence results in Sobolev spaces for a stationary transport equation, Ricerche Mat. Suppl. 36 (1987) 173–184. [8] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [9] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1991. [10] M. Farhloul, Méthodes d’Éléments Finis Mixtes et des Volumes Finis, Thèse, Université Laval, Québec, 1991. [11] M. Farhloul, M. Fortin, A new mixed finite element for the Stokes and elasticity problems, SIAM J. Numer. Anal. 30 (1993) 971–990. [12] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, in: Springer Ser. Comput. Math., Vol. 5, Springer, Berlin, 1986. [13] R.B. Kellogg, B. Liu, A finite element method for the compressible Stokes equations, SIAM J. Numer. Anal. 1 (1996) 780–788. [14] J.R. Kweon, A discontinuous Galerkin method for convection-dominated compressible viscous Navier–Stokes equations with inflow boundary condition, SIAM J. Numer. Anal. 38 (2000) 699–717. [15] J.R. Kweon, An optimal order convergence for a weak formulation of the compressible Stokes system with inflow boundary condition, Numer. Math. 86 (2000) 305–318. [16] J.R. Kweon, A mixed finite element method for a compressible Stokes problem with high Reynolds number, Appl. Numer. Math. 38 (2001) 87–103. [17] J.R. Kweon, R.B. Kellogg, Compressible Navier–Stokes equations in a bounded domain with inflow boundary bondition, SIAM J. Math. Anal. 28 (1997) 94–108. [18] P.D. Lax, R.S. Phillips, Local boundary conditions for dissipative symmetric linear differential operators, Commun. Pure Appl. Math. 18 (1960) 427–455. [19] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, in: Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere, New York, 1980. [20] P.A. Raviart, J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: Lecture Notes in Math., Vol. 606, Springer, Berlin, 1977, pp. 292–315. [21] M.A. Saad, Compressible Fluid Flow, 2nd Edition, Prentice-Hall, London, 1993.
Optimal error estimate for a mixed finite element method for compressible Navier–Stokes system ✩ Jae Ryong Kweon Department of Mathematics, Pohang University of Science and Technology, Pohang 790-784, South Korea
Abstract A linearized stationary compressible viscous Navier–Stokes system is considered. A mixed finite element method is applied and the unique existence of the solution is established by the inf–sup condition. The convection terms, especially in the continuity equation, were thought of causing non-optimal order convergence, but in this paper error estimates of optimal order are derived by implementing the lowest order Raviart–Thomas elements. The error estimates for the normal and tangential components of velocity are also optimal on the interfaces of the interior triangles. It turns out that the non-symmetric discrete system can be reformulated into a symmetric form. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Compressible flows; Mixed finite element; Optimal convergence
1. Introduction So far we have not yet given a complete numerical understanding of the boundary value problems for the Navier–Stokes or the Stokes system for compressible flows [19,21]. Among many unsolved numerical problems is how to construct a finite element solution having optimal error estimate when a finite element method is applied to the system. Note that the system is an elliptic–hyperbolic type. For this, in this paper a mixed finite element method that considers the lowest order Raviart–Thomas elements [8] is applied to a linearized compressible viscous Navier–Stokes equations which can be obtained by linearizing the barotropic stationary compressible viscous Navier–Stokes equations around an ambient flow (U , P ). The equations are −µu − ν∇ div u + ρ (U · ∇)u + ∇p = f in Ω, (1.1) div u + κU · ∇p = g in Ω, u=0 on Γ, ✩
This work was supported by Korea Research Foundation Grant(KRF-2000-015-DP0012). E-mail address: [email protected] (J.R. Kweon).
0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00213-1
276
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
where Ω is a bounded domain with boundary Γ in R2 ; u = [u1 , u2 ] is a velocity vector, p is the pressure; ρ = ρ(P ) is a positive strictly increasing function that provides density as a function of pressure; ρ = dρ/dP , κ = ρ /ρ; the functions U = [U, V ], with U |Γ = 0, and P are the ambient flow which
the functions f and g are given functions, and the numbers µ and ν are the are C 1 functions in Ω; viscous constants with µ > 0 and µ + ν > 0. It is assumed that ρ and κ are given nonzero constants. Otherwise, system (1.1) becomes an incompressible Stokes flow. For the derivation of the equations in (1.1), we refer to [6,13,17]. Note that (1.1) may be obtained by simply adding additional convection terms in the momentum and continuity equations of the (incompressible) Stokes ones. It is often referred to as compressible Navier–Stokes equations. Also note that (1.1) is neither elliptic nor hyperbolic but contains features of each class of equations. In recent years, mixed discretizations have been used in many applications, in particular, for such problems where a major interest lies in the gradient variable rather than the primal one, for example, the flux in stationary incompressible flow problems or neutron diffusion or elasticity problems [2,4,20]. Recently in [16] a mixed scheme was applied to a simple compressible Stokes equation, so an error estimate was shown but not optimal. The convection terms in the equations were thought of causing not optimal order convergence. However, in this paper we show an optimal order error estimate for the solution of the mixed scheme of (1.1) by implementing the lowest order Raviart–Thomas finite elements. In [13] error estimates for the system were given, using an abstract formulation but not optimal. In [15] a weak formulation is proposed for the system with inflow boundary condition and an optimal error estimate is obtained. In [14] a discontinuous Galerkin method was applied to the system. For the stable finite element method, we refer to [3,12,5,9]. It has been observed in [10] that the finite volume methods are very similar to the Raviart–Thomas finite element when they were applied to the Poisson equation. Later in [11] this observation was generalized to the Stokes and elasticity problems by introducing a (new) finite element scheme. This scheme exhibits local properties of conservation of mass and momentum: In fact the method may be viewed as a finite volume scheme. In this paper, based on the integration by part formula (1.3) we formulate a mixed variational scheme and show stability results for the continuous and discrete problems, and also derive optimal error estimates in the variables of the velocity, the gradient of the velocity, and the pressure. Furthermore it is shown that the resulted non-symmetric discrete form for the mixed scheme of (1.1) can be cast into a symmetric matrix form (see (3.62) and (3.64)). Setting the gradient of the velocity u = (u1 , u2 ) by a tensor variable σ , system (1.1) is equivalently written as σ = ∇u in Ω, − div µ σ + ν tr(σ ) − p δ + ρ(U · ∇)u = f in Ω, (1.2) in Ω, tr(σ ) + κU · ∇p = g u=0 on Γ, where σ = (σij )1i,j 2 , σij = ∂ui /∂xj , δ is the identity tensor, and for a tensor τ = (τij )1i,j 2 ,
(div τ )i =
2 ∂τij j =1
∂xj
,
i = 1, 2,
and tr(τ ) = τ11 + τ22 is the trace of the tensor τ .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
277
In formulating a weak formulation of (1.2), the following integration by part formula will be used: σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx = 0. (1.3) Ω
Ω
For showing the stability result, some restrictions are needed. Let µ0 = min{µ, µ + 2ν}, µ1 = µ + 2|ν| and γ∗ = | div U |∞ max{ρ, κ}. In this paper we assume that the numbers µ0 , µ1 and γ∗ satisfy cγ∗ <
µ0 µ21
and
1 C < , µ0 γ∗
(1.4)
where c and C are generic constants. Note that when ν = 0, (1.4) becomes cγ∗ < 1/µ < C/γ∗ . In this paper a finite element scheme for (1.2) is given and unique existences of the continuous and discrete solutions are shown. Our main result is, which is derived in Section 3. Let (σ˜ = (σ, p), u˜ = (u, λn , λt )) and (σ˜ h = (σh , ph ), u˜ h = (uh , λnh , λt h )) be the solutions of (2.11) and (3.1), respectively. Define
κ∗ = Cµ0 µ−2 α∗ = µ(µ + 2ν)/ 3(µ + ν)2 + 5ν 2 . µ∗ = µ0 − Cγ∗ , 1 − γ∗ , If (1.4) holds and µ + 2ν = 0, then the following error estimates hold: ˆ 0,Ω , µ∗ σ − σh 0,Ω + κ∗ p − ph 0,Ω C σ − σˆ 0,Ω + p − p ˆ 0,Ω ), α∗ u − uh 0,Ω C(σ − σh 0,Ω + u − u
(1.5)
where (σˆ , p) ˆ ∈ Xh and uˆ ∈ V 1h are the interpolations of (σ, p) and u, defined by (3.9) and (3.12), respectively. Note that the numbers µ∗ and κ∗ are well-defined under (1.4), and also α∗ = 0 with µ + 2ν = 0. Note that (1.5) is optimal in a sense that the errors in both sides were measured with the same norms. We give the error estimate for the normal and tangential components of velocity on the interelement boundaries. Let λn = u · n and λt = u · t be the normal and tangential components of u on the boundary ∂K of the triangle K ∈ Th where Th is a given triangulation. Assume that µ + 2ν = 0 and (1.4) holds. If div U is constant on each triangle K ∈ Th , then 0 ρ (u · t) − λt h + ρh0 (u · n) − λnh 1/q,p,∂K Ch2/p Z1 , (1.6) h 1/q,p,∂K where Z1 = |σ |1,Ω + |p|1,Ω + | div(µ σ + ν tr(σ )δ)|1,Ω and ρh0 is the projection operator defined in (3.14). For k 0 and 1 p < ∞ we denote by W k,p (Ω) the Sobolev space with norm wk,p,Ω and seminorm |w|k,p,Ω (see [1]). Let W 1/q,p (Γ ) be the space of the traces v|Γ over Γ of the functions 1,p v ∈ W 1,p (Ω) where 1/p + 1/q = 1. Also W0 (Ω) denote the space of functions v in W 1,p with 1,p zero boundary value. W −1,q (Ω) and W −1/q,q (Γ ) denote the dual spaces of W0 (Ω) and W 1/q,p (Γ ), 1,2 respectively. If p = 2, W k,2 (Ω) and W0 (Ω) are denoted by H k (Ω) and H01 (Ω), respectively, and the norm uk,2,Ω is replaced by uk,Ω . Also L20 (Ω) denote the space of functions in L2 (Ω) with zero mean value over Ω. We denote by
2 Wq (div; Ω) = v ∈ Lq (Ω) ; div v ∈ Lq (Ω) , where 1 q < ∞, and W2 (div; Ω) means H(div; Ω). Finally we set
Mβ = χ ∈ L2 (Ω): χβ,Ω < ∞ ,
278
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
with norm χβ,Ω = χ20,Ω + div(U χ)20,Ω . We use the Sobolev embedding theorems: W 1,p (D) 1→ Lq (D) (p q 2p/(2 − p)), and W 1/q,p (D) 1→ L2 (∂D) (4/3 < p < q), where D ⊂ R 2 is a bounded domain. For the tensor τ = (τij )1i,j 2 we denote by τi , i = 1, 2, the vector τi = (τi1 , τi2 ) and the normal trace τ n = (τ1 · n, τ2 · n). The normal and tangential components, τnn and τnt , of the vector τ n shall be defined by τnn = τ n · n =
2
τ1j nj n1 +
j =1
τnt = τ n · t =
2 j =1
τ1j nj t1 +
2
τ2j nj n2 ,
j =1 2
τ2j nj t2 ,
j =1
where t denotes the unit tangential vector to Γ. The paper is organized as follows. In Section 2, based on the integration by part formula, several bilinear forms are defined and a variational formulation for (1.2) is presented, which consider the Lagrange multipliers on inter-element boundaries and an inf–sup condition is shown for the stability of the continuous problem. In Section 3 the continuous problem is discretized and the stability result is shown by the discrete inf–sup condition. Finally optimal error estimates are derived and the nonsymmetric matrix generated by the discrete problem is reformulated into a symmetric form. We use C to denote a generic positive constant, not depending on the mesh size h. Note that C may take different values in different places.
2. Mixed formulation In this section we first give a weak form for (1.2) (see (2.2)), where in [16] a similar form was analyzed, and consider the weak form that employs the Lagrange multipliers on the interelement boundaries (see (2.6)) and show a unique existence. We next give a regularity result for (1.1), which can be derived from [6]: Theorem A. Let Ω be a bounded domain of class C 2 and q > 2. Assume that f ∈ (Lq (Ω))2 and g ∈ W 1,q (Ω) with Ω g dx = 0. Then there is a constant K1 > 0 such that if U 2,q,Ω < K2 for a constant K2 , there is a unique solution (u, p) ∈ (W 1,q (Ω))2 × W 0,q (Ω) of (1.1), satisfying u2,q,Ω + p1,q,Ω K1 (f 0,q,Ω + g1,q,Ω ). For the weak formulation of (1.2), we define Qβ = Mβ ∩ L20 (Ω), V β = (Mβ )2 and
2 X0 = (τ, χ) ∈ (V 0 )2 × Q0 : µτ + ν tr(τ ) − χ δ ∈ H (div; Ω) .
(2.1)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
279
Note that V0 = L2 (Ω) and Q0 = L20 (Ω). Using (1.3) and (2.1), a weak formulation of (1.2) reads: find ((σ, p), u) ∈ X0 × V 0 such that σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx = 0, Ω Ω (2.2) div µσ + ν tr(σ ) − p δ · v + ρu · div(U v) dx = − f · v dx, Ω Ω tr(σ )χ − κp div(U χ) dx = gχ dx, Ω
Ω
for all (τ, χ) ∈ X0 , v ∈ V β and χ ∈ Qβ . It is obvious that (2.2) is equivalent to (1.2). In this paper we consider a finite element discretization for (2.2) by introducing the Lagrange multipliers on the interelement boundaries so that the normal trace (µτ + (ν tr(τ ) − χ)δ)n is continuous
This Lagrange multiplier is nothing else than the across interelement edges of the triangulation Th of Ω. normal component u · n and the tangential component u · t of the velocity. The tensor variable will be approximated by the lowest degree of the Raviart–Thomas finite element [8] 2 RT0 (K) := P0 (K) + x · P0 (K), where P0 (K) is the space consisting of constant functions on K and x = (x1 , x2 ) ∈ R2 , while the velocity and the pressure will be constant on each element of the triangulation and the Lagrange multiplier will be constant on every edge of each element. Note that the space (H (div; K))2 , K ∈ Th , is not appropriate in obtaining an approximation of the tensor τ by the Raviart–Thomas finite elements. In fact, if an interpolation operator is to be defined on the space (RT0 (K))2 by using the three normal components of a tensor τ on ∂K, K ∈ Th as the degrees of freedom, then the tensor τ should be slightly smoother than merely belonging to (H (div; K))2 [8]. From this point of view we will require q > 2 in the index q of the space (Wq (div; K))2 , which can be considered a finite element space for the tensor variable. Now we will choose the corresponding functional spaces that the Lagrange multipliers can be defined on interelement edges of the triangulation Th and the normal trace (µτ + (ν tr(τ ) − χ)δ)n is continuous there. Let p be a real number with 4/3 < p < 2 and q its conjugate where 1/p + 1/q = 1. We consider the following spaces [11]: Let V 1 = (Lq (Ω))2 , M = L20 (Ω) ∩ Lq (Ω) and
2 X = (τ, χ) ∈ V 1 × M: µτ + ν tr(τ ) − χ δ |K ∈Wq (div; K), ∀K ∈ Th , (2.3) = η ∈ W 1/q,p (∂K): η|∂K∩Γ = 0, ∀K ∈ Th . V 2 K∈Th
1 = V β ∩ (Lq (Ω))2 and We also denote by V V = V 1 × (V2 )2 ,
=V 1 × (V2 )2 , V
Q = Qβ ∩ Lq (Ω),
(2.4)
where Vβ and Qβ are defined above. The following variables are introduced: σ˜ = (σ, p),
τ˜ = (τ, χ),
u˜ = (u, λn , λt ),
v˜ = (v, ηn , ηt ).
(2.5)
280
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using the spaces in (2.3), (2.4) and the variables given in (2.5), a weak formulation for (1.2) reads: find σ˜ = (σ, p) ∈ X, u˜ = (u, λn , λt ) ∈ V such that σ : µτ + ν tr(τ ) − χ δ dx + div µτ + ν tr(τ ) − χ δ · u dx Ω K∈Th K τnt λt ds = 0, ∀τ˜ ∈ X, µτ + ν tr(τ ) − χ δ nn λn ds − µ − ∂K ∂K K∈ T K∈ T h h 1, div µσ + ν tr(σ ) − p δ · v + ρu div(U v) dx = − f v dx, ∀v ∈ V K Ω K∈Th (2.6) η ds = 0, ∀η ∈ V , µσ + ν tr(σ ) − p δ n n 2 nn K∈Th ∂K σnt ηt ds = 0, ∀ηt ∈ V2 , µ ∂K K∈ T h χ tr(σ ) − κp div(U χ) dx = gχ dx, ∀χ ∈ Q. K Ω K∈Th
Note that the Lagrange multipliers were introduced in (2.6) but not in (2.2). In order to write (2.6) in a compact form we consider four bilinear forms on the functional spaces in (2.3) and (2.4): (2.7) a(σ˜ , τ˜ ) = σ : µτ + ν tr(τ ) − χ δ dx, σ˜ , τ˜ ∈ X, Ω
and b(τ˜ , v˜ ) =
K∈Th K
−µ
div µτ + ν tr(τ ) − χ δ · v dx − µτ + ν tr(τ ) − χ δ nn ηn ds
K∈Th∂K
τnt ηt ds,
τ˜ ∈ X, v˜ ∈ V ,
(2.8)
K∈Th∂K
and c(p, χ) = −
p div(U χ) dx,
p ∈ M, χ ∈ Q,
(2.9)
K∈Th K
and d(σ, χ) =
χ tr(σ ) dx,
σ ∈ (V 1 )2 , χ ∈ M.
(2.10)
K∈Th K
Using (2.7)–(2.10), (2.6) can be written by: find σ˜ ∈ X and u˜ ∈ V such that ˜ = 0, τ˜ ∈ X, a(σ˜ , τ˜ ) + b(τ˜ , u) , b(σ˜ , v˜ ) − c(ρ u, v) = −f , v , ∀v˜ ∈ V c(κp, χ) + d(σ, χ) = g, χ , ∀χ ∈ Q, where u, v = Ω uv dx.
(2.11)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
281
Note that the form a(σ˜ , τ˜ ) is not symmetric on X × X. However, this non-symmetric form a can be decomposed into a symmetric form plus the other bilinear one by taking the test function τ˜ ∈ X suitably. For details, see (3.58). In order to show the unique existence of solution of (2.11), which is shown in Theorem 2.1 below, we recall the inf–sup condition for the bilinear form b, which can be shown by the same argument as given in [16, Lemma 2.1]: There is a positive constant α such that sup τ˜ ∈X τ˜ =0
b(τ˜ , v˜ ) α v0,Ω , τ˜ X
∀v˜ ∈ V ,
(2.12)
where τ˜ X := (τ 20,Ω + χ20,Ω )1/2 . Now using (2.12) and Theorem A, we show a unique existence of the solution for the continuous problem (2.11): Theorem 2.1. Assume that the conditions given in Theorem A hold. Assume that β2,q,Ω is small in a suitable sense. Let (u, p) be the solution of (1.1) and σ = ∇u. Then (σ˜ = (σ, p), u˜ = (u, u · n, u · t)) is the unique solution of (2.11). Furthermore, if (1.4) holds, then µ∗ σ 0,Ω + κpβ 0,Ω + κ∗ p0,Ω C f 0,Ω + g0,Ω , √ (2.13) αu0,Ω 2(µ1 + 1)σ 0,Ω , √ where µ∗ = µ0 − Cαγ∗ and κ∗ = Cµ0 µ−2 2 − γ∗ with µ2 = µ1 + ρ|U |∞ . Proof. Since the conditions given in Theorem A is assumed, the pair σ˜ = (σ, p) ∈ X, u˜ = (u, u · n, u · t) ∈ V and is the solution of (2.11). The uniqueness follows from Theorem A. The inequality (2.13), constants µ∗ and κ∗ can be derived by applying the same procedures as given in the proof of [16, Theorem 2.1]. ✷
3. Mixed discretization Before proceeding further, we recall that in [16] a mixed finite element scheme that implements the Raviart–Thomas elements of order l was applied to a (simple) compressible Stokes system but does not consider the Lagrange multipliers on the interelement boundaries like (2.6). However, since the governing equations and finite elements considered in this paper have generically similar characters like [16], we often use some similar arguments as employed in [16] due to technicality and complexity. Now we are going to approximate the solution of (2.11) by implementing the lowest order Raviart– Thomas finite elements. Let h be an approximate indicator. With finite-dimensional spaces Xh ⊂ X, V h ⊂ V , Mh ⊂ M and Qh ⊂ Q, we find the approximate solutions σ˜ h ∈ X h and u˜ h ∈ V h such that ∀τ˜h ∈ X h , a(σ˜ h , τ˜h ) + b(τ˜h , u˜ h ) = 0, (3.1) b(σ˜ h , v˜ h ) − c(ρ uh , vh ) = −f , vh , ∀v˜ h ∈ V h , ∀χh ∈ Qh . c(κph , χh ) + d(σ˜ h , χh ) = g, χh , In order to analyze the discrete problem (3.1) we are going to specify the finite-dimensional spaces Xh , V h , Mh , Qh and their degrees of freedom. Let Th be a given regular triangulation for the closure of Ω and Pl (K) the space of polynomials of degree l, where K ∈ Th is a triangle.
282
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Let RT0 (K) denote by the lowest order Raviart–Thomas element 2 RT0 (K) := P0 (K) + x · P0 (K). Note that any v ∈ RT0 (K) is uniquely determined by its normal components n · v |ei on the edges ei , 1 i 3, of K ∈ Th , where n denotes the outer normal vectors of K. Letting D := (RT0 (K))2 , we define an operator ΠK1 : (Wq (div; K))2 → D by ΠK1 τ ∈ (RT0 (K))2 such that 2 1 (3.2) ΠK τ − τ n · v 0 ds = 0, ∀v 0 ∈ R0 (∂K) , ∂K
where R0 (∂K) = {v: v|e ∈ P0 (e), ∀e ∈ ∂K}. Note that the boundary integral in (3.2) may be written by 1 (3.3) ΠK τ n · v 0 ds = τ n, v0 ∂K , ∂K
where . , . ∂K denotes the duality pairing between (W −1/q,q (∂K))2 and (W 1/q,p (∂K))2 . The interpolation error estimates for ΠK1 are cited in [11]: if Th is regular, then 1 Π τ − τ ChK |τ |1,K , K 0,K 1 div Π τ − τ C div τ 0,K , K
0,K
(3.4) (3.5)
where hK is a diameter of K. Using the operator ΠK1 of (3.3), let Πh1 be defined by Πh1 τ|K = ΠK1 τ,
∀K ∈ Th ,
(3.6)
for all τ ∈ {τ ∈ (L (Ω)) : τ ∈ (Wq (div; K)) , ∀K ∈ Th }. Define an approximate space of X as follows:
(3.7) X h = τ˜h ∈ X: τh|K ∈ D, χh|K ∈ P0 (K), ∀K ∈ Th , q
4
2
where τ˜h = (τh , χh ). Let ρh be an projection operator defined by 1 rh χ dx, ρh χ = rh χ − |Ω|
(3.8)
Ω
where rh is the L2 projection onto the space consisting of piecewise constants. Noting that (µ τh + (ν tr(τh ) − χh )δ)|K ∈ (RT0 (K))2 for ∀K ∈ Th , we define an operator Πh : X → X h , Πh τ˜ = (τh , χh ) ∈ X h (see [11]) such that χh = ρh χ and for ∀v ∈ (R0 (∂K))2 and ∀K ∈ Th ,
µτ + ν tr(τ ) − χ δ − µτh + ν tr(τh ) − χh δ n · v ds = 0. (3.9) ∂K
Using (3.7) and (3.9) we have µτh + ν tr(τh ) − χh δ = Πh1 µτ + ν tr(τ ) − χ δ , and the interpolation error estimate for Πh is given in [11]: assuming that Th is regular and (τ, χ) ∈ (H 1 (Ω))4 × H 1 (Ω), we have (3.10) Πh τ˜ − τ˜ 0,Ω Ch |τ |1,Ω + |χ|1,Ω .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
For the approximations of V , M and Q, we set
2 V 1h = v h ∈ V 1 : v h |K ∈ P0 (K) , ∀K ∈ Th ,
V2h = ηh ∈ V2 : ηh |∂K ∈ R0 (∂K), ∀K ∈ Th , V h = V 1h × (V2h )2 , M = χ ∈ M: χ | ∈ P (K), ∀K ∈ T }, Qh = Mh ∩ Q. h h K 0 h Define the operator Ph1 : V 1 → V 1h by 2 Ph1 v|K = PK1 v, ∀v ∈ V 1 , PK1 v ∈ P0 (K) such that 2 (PK1 v − v) · v 0 dx = 0, ∀v 0 ∈ P0 (K) , ∀K ∈ Th .
283
(3.11)
(3.12)
K
Consider the operator Ph : V → V h defined by Ph v˜ = Ph (v, ηn , ηt ) = Ph1 v, ρh0 ηn , ρh0 ηt , where ρh0 is the projection operator from K∈Th L2 (∂K) onto K∈Th R0 (∂K), i.e., 0 0 ρh : V2 → V2h , ρh η ∈ R0 (∂K), 0 ρh η − η v0 ds = 0, ∀v0 ∈ R0 (∂K), ∀K ∈ Th .
(3.13)
(3.14)
∂K
Note that from the Sobolev’s imbedding theorem W 1/q,p (∂K) 1→ L2 (∂K)(4/3 < p < q), the operator ρh0 is well-defined on V2 (cf. see [11]). In next lemma we establish a discrete inf–sup condition for the bilinear form b, in which the proof is similar to the ones in [11] and [16, Lemma 3.1]: Lemma 3.1. Assume that Th is regular and µ + 2ν = 0. Then there exists a constant α∗ independent of h such that b(τ˜h , v˜ h ) α∗ v h 0,Ω , ∀v˜ h ∈ V h , (3.15) sup τ˜h ∈X h τ˜h X −1/2
where α∗ = C/µ3
with µ3 defined in (3.16).
Proof. Let v˜ h = (v h , µnh , µt h ) ∈ V h . By [11, Proposition 3.1] there is a τh∗ ∈ (Wq (div; Ω))2 such that τh∗ |K ∈ (RT0 (K))2 , ∀K ∈ Th , div τh∗ = v h in Ω and τh∗ 0,Ω Cv h 0,Ω . Picking a function τh satisfying µτh + ν tr(τh )δ = τh∗ , the components of τh can be expressed in terms of those of τh∗ and (µ + ν)2 + ν 2 ∗ 2 ∗ 2 + τ τ h11 0,Ω h22 0,Ω (µ2 + 2µν)2 1 ∗ 2 ∗ 0,Ω + τh21 20,Ω µ3 τh∗ 20,Ω Cµ3 v h 20,Ω , + 2 τh12 µ
τh 20,Ω 2
(3.16)
where µ3 = (3(µ + ν)2 + 5ν 2 )/(µ2 (µ + 2ν)2 ). Finally, letting τ˜h = (τh , 0) ∈ X h , we have b(τ˜h , v˜ h ) = −1/2 v h 20,Ω (since τh ∈ (Wq (div; Ω))2 ). Thus the required inequality (3.15) follows with α∗ = Cµ3 . ✷
284
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using Lemma 3.1 and considering some similar arguments used in the proof of [16, Theorem 3.1], we show a unique existence of the discrete solution of (3.1). Theorem 3.1. Let µ + 2ν = 0. Assume that (1.4) holds. Then problem (3.1) has a unique solution (σ˜ h , u˜ h ) ∈ X h × V h , satisfying µ∗ σh 0,Ω + κ∗ ph 0,Ω C f 0,Ω + g0,Ω , (3.17) α∗ uh 0,Ω Cσh 0,Ω , where µ∗ and κ∗ are given in (3.26), and α∗ is given in (3.15), and C is a constant not depending on h. Proof. Let Zh = {τ˜h ∈ X h : b(τ˜h , v˜ h ) − c(ρuh , v h ) = −f , vh , ∀v˜ h ∈ V h } and let τ˜h = (τh , χh ) ∈ Zh be arbitrary. Recalling that V 1h consists of piecewise constant functions, and letting (div U )h and f h be the L2 projections into the piecewise constant functions which are constants on each triangle K ∈ Th , we have (3.18) div µτh + ν tr(τh ) − χh δ + ρuh (div U )h + f h v h dx = 0, K∈Th K
for all v˜ h = (vh , 0, 0) ∈ V h . It follows from (3.18) that div µτh + ν tr(τh ) − χh δ + ρuh (div U )h = −f h .
(3.19)
Setting σ1 = µτh + (ν tr(τh ) − χh )δ with χh ∈ L20 (Ω), we have div σ1 = −ρuh (div U )h − f h . We now need to show that (3.20) χh 0,Ω c1 µ1 τh 0,Ω + c2 γ∗ uh 0,Ω + f h 0,Ω , ∀τ˜h ∈ Zh , √ √ and c2 = C/ 2. Indeed, to show (3.20), we define where c1 = 2(2 + C) µ + 2ν |Ω|−1 tr(τh ) dx, L= 2 Ω
and σ2 = µτh + (ν tr(τh ) − χh − L)δ. Using (3.19), we have div σ2 = −ρuh (div U )h − f h , and Ω tr(σ2 ) dx = 0. But using [8, p. 161, Proposition 3.1], there is a constant C such that σ2D + div σ2 0,Ω , (3.21) σ2 0,Ω C 0,Ω where σ2D = σ2 − 12 tr(σ2 )δ. On the other hand, µ σ2D = µτh − tr(τh )δ, 2 and using the similar procedures as given in [16, Theorem 3.1], h 0,Ω + C div σ2 0,Ω + Lδ0,Ω σ1 0,Ω µCτ + 1 + 2|ν| τh 0,Ω + C γ∗ uh 0,Ω + f h 0,Ω . µ 2C Since χh δ = µτh + ν tr(τh )δ − σ1 , and using (3.23) we have 1 χh 0,Ω √ µ1 τh 0,Ω + σ1 0,Ω c1 µ1 τh 0,Ω + c2 γ∗ uh 0,Ω + f h 0,Ω , 2
(3.22)
(3.23)
(3.24)
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
285
√ √ and c2 = C/ 2. Thus (3.20) follows. Now, recalling µ0 = min{µ, µ + 2ν} and where c1 = 2(2 + C) using (3.20), we obtain f , uh + g, ph = a(σ˜ h , σ˜ h ) + d(σh , ph ) + c(ρuh , uh ) + c(κph , ph ) µ0 σh 20,Ω − γ∗ uh 20,Ω + ph 20,Ω µ0 µ0 −1 2 − α∗ γ∗ σh 0,Ω + − γ ph 20,Ω 2 4(c1 µ1 )2 µ0 2 2 2 − u + f γ h h ∗ 0,Ω 0,Ω 2(c1 µ1 )2 1 (3.25) µ∗ σh 20,Ω + κ∗ ph 20,Ω − Cf h 20,Ω 2 where µ∗ = (µ0 − c3 γ∗ ) with c3 = 2α∗−1 + µ0 γ∗ /(c1 µ1 )2 , and κ∗ = µ0 /4(c1 µ1 )2 − γ∗ . Finally, computing f , uh + g, ph and using f h 0,Ω f 0,Ω , the inequality (3.17a) follows from (3.25). Also (3.17b) follows from Lemma 3.1 and Eq. (3.1a). The uniqueness follows from (3.17). ✷ From the definition of the operators Πh and Ph by (3.9) and (3.13), respectively, we recall that the bilinear form b defined by (2.8) satisfies b(Πh τ˜ − τ˜ , v˜ h ) = 0, b(τ˜h , Ph v˜ − v˜ ) = 0,
∀v˜ h ∈ V h , ∀τ˜ ∈ X, ∀v˜ ∈ V , ∀τ˜h ∈ X h .
(3.26)
We are now ready to derive an optimal error estimate which can be compared with [16, Theorems 3.2 and 4.2]. The derivation is similar to the one in [16, Theorems 3.2]. Theorem 3.2. Assume that (1.4) holds and µ + 2ν = 0. Let Th is a given regular triangulation of the domain Ω. Let (σ˜ , u) ˜ be the solution of (2.11) and (σ˜ h , u˜ h ) the solution of (3.1). Assume that 1 4 σ˜ = (σ, p) ∈ (H (Ω)) × (H 1 (Ω) ∩ L20 (Ω)). Then there exists a constant C independent of h such that ˆ 0,Ω + u − u ˆ 0,Ω , µ∗ σ − σh 0,Ω + κ∗ p − ph 0,Ω C σ − σˆ 0,Ω + p − p (3.27) ˆ 0,Ω , α∗ u − uh 0,Ω C σ − σh 0,Ω + u − u where σˆ , pˆ and uˆ are the interpolations of σ, p and u, respectively. Proof. From (2.11) and (3.1), and using (3.26), we have ˜ = a(σ˜ − Πh σ˜ , τ˜h ), ∀τ˜h ∈ X h , a(σ˜ h − Πh σ˜ , τ˜h ) + b(τ˜h , u˜ h − Ph u) 1 b(σ˜ h − Πh σ˜ , v˜ h ) − c uh − Ph u, vh = −c u − Ph1 u, v h , ∀v˜ h ∈ V h , c(ph − ρh p, χh ) + d σh − Πh1 σ, χh = c(p − ρh p, χh ) + d σ − Πh1 σ, χh , Applying to (3.29) the same procedures used in showing (3.20), we have ˆ 0,Ω c1 µ1 σh − Π 1 σ 0,Ω + c2 γ∗ uh − P 1 u + u − P 1 u ph − p h
h
0,Ω
h
(3.28) (3.29) ∀χh ∈ Qh .
, 0,Ω
(3.30)
(3.31)
where ci (i = 1, 2) are given in (3.20). Taking v˜ h = u˜ h − Ph u˜ in (3.28), τ˜h = σ˜ h − Πh σ˜ in (3.29), and χh = ph − pˆ with pˆ = ρh p in (3.30), and adding the resulted equations, and applying the same procedures used in showing (3.25), and using (3.31), we obtain
286
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
a(σ˜ h − Πh σ˜ , σ˜ h − Πh σ˜ ) + d σh − Πh1 σ, ph − pˆ ˆ ph − p) ˆ + c uh − Ph1 u, uh − Ph1 u + c(ph − p, 2 2 1 ˆ 20,Ω − C u − Ph1 u 0,Ω , (3.32) µ∗ σh − Πh1 σ 0,Ω + κ∗ ph − p 2 where µ∗ and κ∗ are defined in (3.25). On the other hand, the resulting right-hand is estimated by a(σ˜ − Πh σ˜ , σ˜ h − Πh σ˜ ) + c u − Ph1 u, uh − Ph1 u ˆ + d σ − Πh1 σ, ph − pˆ + c(p − p, ˆ ph − p) 2 2 ˆ 20,Ω + ε3 uh − Ph1 u 0,Ω ε1 σh − Πh1 σ 0,Ω + ε2 ph − p 2 2 + C σ − Πh1 σ 0,Ω + p − p ˆ 20,Ω + u − Ph1 u 0,Ω , ∀εi > 0, i = 1, 3, (3.33) where C = C(γ∗ , µ1 , |U|∞ ). Furthermore, using Lemma 3.1 and (3.28), , α∗ uh − P 1 u C σh − Π 1 σ + u − P 1 u h
h
0,Ω
0,Ω
h
0,Ω
(3.34)
where α∗ is given in (3.15). Thus, combining (3.32)–(3.34) and using the triangle inequality, the required inequality follows. ✷ In the formulation of (2.11), λn (respectively, λt ) is nothing else than the normal (respectively, the tangential) component u · n (respectively u · t) of the velocity. Therefore, in the discrete problem (3.1), λnh (respectively, λt h ) is an approximation of u · n (respectively, u · t). We are now going to derive error bounds for λnh − u · n and λt h − u · t. To do this, we cite the following lemma (see [11]): Lemma 3.2. For any T ∈ (W −1/q,q (∂K))2 , there exists a unique τh ∈ (RT0 (K))2 such that 2 (T − τh n) · v0 ds = 0, ∀v0 ∈ R0 (∂K) , ∂K
where R0 (∂K) = {v: v|e ∈ P0 (e), ∀e ∈ ∂K}. Furthermore, assume that Th is uniformly regular. Then there is a constant C independent of h such that τh 0,q,K CT −1/q,q,∂K , τh 0,K Ch2/p−1T −1/q,q,∂K , div τh 0,K Ch2/p−2T −1/q,q,∂K , where p and q are conjugate. Using Lemma 3.2 and applying exactly the same arguments given in [11, Theorem 3.3], we obtain: Theorem 3.3. Let (σ, u) be the solution of (2.11) and (σ˜ h , u˜ h ) the solution of (3.1). Assume that µ + 2ν = 0. If Th is uniformly regular and σ˜ = (σ, p) ∈ (H 1 (Ω))4 × (H 1 (Ω) ∩ L20 (Ω)), then there exists a constant C independent of h such that 0 ρ (u · t) − λt h + ρ 0 (u · n) − λnh Ch2/p−1 Z, (3.35) h
1/q,p,∂K
h
for all K ∈ Th , where Z = |σ |1,Ω + |p|1,Ω .
1/q,p,∂K
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
287
The error estimate (3.35) is not optimal. In order to get an optimal estimate we need the following Lemmas: Lemma 3.3. (1) Assume that |∇ U |∞ is small enough. Then there is a unique solution z of the problem: µz + ν∇ div z + ρ div(U z) = w z=0
in Ω, on Γ,
(3.36)
satisfying µz1,Ω Cw0,Ω , with C = C(|∇U |∞ ). Also the solution z satisfies µzk,q,Ω Cwk−2,q,Ω with C = C(Ω, U , k), k is an integer 1, 1 < q < ∞. (2) Let k −1 be an integer and let q ∈ (2/(k + 2), ∞). Assume that U ∈ W k+3,q (Ω) with U |Γ = 0 and Γ ∈ C k+3 . Assume that µ + 2ν = 0 and 2/(µ + 2ν) > αk ≡ CU k+3,q,Ω . For any S ∈ W k,q (Ω) there is a solution s ∈ W k,q (Ω) of problem 2 s − κ div(U s) = S, µ + 2ν
(3.37)
2 − αk )sk,q,Ω CSk,q,Ω for a constant C. satisfying ( µ+2ν
Proof. (1) It easily follows from the standard Lq theory of elliptic equations and [17, Lemma 2.4]. (2) It follows from [7, Theorem 1.1] that problem (3.37) has a solution s satisfying the required inequality. For further references one may refer to Beirão da Veiga [6, p. 234] and Lax & Phillips [18]. ✷ On the basis of (3.36) and (3.37) we consider the following problem µz + ν∇ div z + ρ div(U z) = w in Ω, 2 s − κ div(U s) + div z = h in Ω, µ+2ν z=0 on Γ, s dx = 0,
(3.38)
Ω
where z and s are the unknown variables, and w, h are given functions. From Lemma 3.3 one can derive that if w ∈ L2 (Ω) and h ∈ H 1 (Ω) are given, then there is a unique solution (z, s) ∈ (H01 (Ω))2 × L20 (Ω) of (3.38) satisfying z2,Ω K1 w0,Ω , s1,Ω K2 h1,Ω + w0,Ω , where K1 is a constant and K2 = C(2/(µ + 2ν) − α1 )−1 . In order to derive a dual formulation of (1.2), let γ = (γij )1i,j 2 be a tensor variable and let µγ + ν tr(γ ) − s δ = µ∇z + ν div zδ,
(3.39) (3.40)
(3.41)
where δ is the identity tensor. Taking the trace to both sides of (3.41), we have tr(γ ) =
2 s + div z. µ + 2ν
(3.42)
288
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Using (3.41) and (3.42), system (3.38) becomes µγ + ν tr(γ ) − s δ = µ∇z + ν div zδ div µγ + ν tr(γ ) − s δ + ρ div(U z) = w tr(γ ) − κ div(U s) = h z=0 s dx = 0.
in Ω, in Ω, in Ω, on Γ,
(3.43)
Ω
On the other hand, multiplying τ to both sides of (3.41) and doing the integration by part and using z|Γ = 0, we have the following formula (see (1.3)): τ : µγ + ν tr(γ ) − s δ dx + div µτ + ν tr(τ )δ · z dx = 0. Ω
Ω
Hence a dual formulation for (2.11) is: Find γ˜ = (γ , s) ∈ X, z˜ = (z, αn , αt ) ∈ V such that τ˜ ∈ X, a(τ˜ , γ˜ ) + b1 (τ, z˜ ) = 0, b(γ˜ , v˜ ) − c(ρv, z) = (w, v), ∀v˜ ∈ V , c(κ χ, s) + d(γ , χ) = (h, χ), ∀χ ∈ M, where div µτ + ν tr(τ )δ · z dx − µτ + ν tr(τ )δ nn αn + µτnt αt ds b1 (τ, z˜ ) = K∈Th K
(3.44)
K∈Th∂K
for τ˜ = (τ, χ) ∈ X and z˜ ∈ V ,
(3.45)
and the bilinear forms a, b, c and d are defined in (2.7)–(2.10), respectively. Using (3.39), (3.40) and (3.43c) we obtain the following lemma: Lemma 3.4. Let µ + 2ν = 0. Assume that w ∈ (L2 (Ω))2 and h ∈ H 1 (Ω). Assume that the solution (z, s) of (3.38) has the regularity µγ + (ν tr(γ ) − s)δ ∈ (H (div; Ω))2 and s ∈ L20 (Ω), where µγ + (ν tr(γ ) − s)δ = µ∇z + ν div zδ. Then γ˜ = (γ , s) and z˜ = (z, z · n, z · t) is the unique solution of (3.44) and also satisfy z2,Ω Cw0,Ω , γ 1,Ω + s1,Ω C h1,Ω + w0,Ω ,
(3.46) (3.47)
where C = C(K1 , K2 ). Using the estimates (3.46) and (3.47), optimal error estimates for u · n and u · t are derived: Theorem 3.4. Assume that the conditions of Theorem 3.3 hold and div[µσ + ν tr(σ )δ] ∈ (H 1 (Ω))2 . Assume that µ + 2ν = 0 and (1.4) hold. Then if div U is constant on each triangle K ∈ Th , there exists a constant C not depending on h such that 0 ρ (u · t) − λt h + ρ 0 (u · n) − λnh Ch2/p Z1 , h
1/q,p,∂K
h
1/q,p,∂K
where Z1 = |σ |1,Ω + |p|1,Ω + | div[µσ + ν tr(σ )δ]|1,Ω .
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
Proof. Let γ˜ = (γ , s) ∈ X, z˜ = (z, αn , αt ) ∈ V be the solution of ∀τ˜ = (τ, χ) ∈ X, a(τ˜ , γ˜ ) + b1 (τ, z˜ ) = 0, b(γ˜ , v˜ ) − c(ρ v, z) = (u¯ h − uh , v), ∀v˜ ∈ V , c(κ χ, s) + (tr(γ ), χ) = 0, ∀χ ∈ M
289
(3.48)
where u¯ h = Ph1 u (Ph1 is defined by (3.12)). Using u¯ h − uh ∈ (L2 (Ω))2 and Lemma 3.4, we have z1,Ω Cu¯ h − uh 0,Ω
and
γ˜ 1,Ω Cu¯ h − uh 0,Ω ,
(3.49)
where C is a positive constant not depending on h. First, noting that c(u¯ h − uh , z) = 0 (by the integration by parts and U |Γ = 0), using (3.48b) and (3.26), we have u¯ h − uh 20,Ω = b(Πh γ˜ , u˜ − u˜ h )
(using (3.28))
= a(σ˜ h − σ˜ , Πh γ˜ ) = a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + a(σ˜ h − σ˜ , γ˜ )
(using (3.48a))
= a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + b1 (σ − σh , z˜ ).
(3.50)
On the other hand, letting Ph z˜ = (¯zh , αhn , αht ), where z¯ h := Ph1 z and αh := z¯ h|∂K , we have b(Πh σ˜ − σ˜ h , Ph z˜ ) = b1 Πh1 σ − σh , Ph z˜ = b1 (σ − σh , Ph z˜ ). Hence it follows from (3.29) and (3.26) that b1 (σ − σh , Ph z˜ ) − ρc u − uh , Ph1 z = 0.
(3.51)
If div U is constant on each triangle K ∈ Th , then c(u − uh , Ph1 z) = c(u¯ h − uh , Ph1 z). Subtracting (3.51) from (3.50), we have u¯ h − uh 20,Ω = a(σ˜ h − σ˜ , Πh γ˜ − γ˜ ) + b1 (σ − σh , z˜ − Ph z˜ ) + ρ c(u¯ h − uh , Ph1 z) (using (3.49)) Ch2 Z1 |γ˜ |1,Ω + |z|1,Ω + C| div U |∞ u¯ h − uh 20,Ω ,
(3.52)
where Z1 = |σ |1,Ω + |p|1,Ω + | div(µ σ + ν tr(σ )δ)|1,Ω and for (3.52) we used div µσ + ν tr(σ )δ − Πh1 µσ + ν tr(σ )δ · z − Ph1 z dx, b1 (σ − σh , z˜ − Ph z˜ ) = K∈Th K
which follows from the properties of the operator Ph and the fact µσ + ν tr(σ )δ ∈ (Wq (div; Ω))2 . Using (3.52) and (3.49) we have u¯ h − uh 0,Ω C1 h2 Z1 ,
(3.53)
where C1 = C/(1 − C| div U |∞ ). Using exactly the same procedures as ones given in [11, Theorem 3.4] and the inequality (3.38) given in [11, Theorem 3.3], we obtain the required inequality. ✷
290
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
As mentioned in Section 2 the form a(σ˜ , τ˜ ) is not symmetric. However, by taking the test function τ˜h ∈ X suitably, system (2.11) can be transformed into a symmetric system. For τ˜ = (τ, 0) ∈ X let ˜ = b(τ˜ , u) ˜ and for τ˜ = (0, χ) ∈ X let a2 (σ˜ , τ˜ ) = a(σ˜ , τ˜ ), b2 (τ˜ , u) ˜ = b(τ˜ , u). ˜ a1 (σ˜ , τ˜ ) = a(σ˜ , τ˜ ), b1 (τ˜ , u) Then a1 (σ˜ , τ˜ ) is a symmetric form for τ˜ = (τ, 0) ∈ X and a2 (σ˜ , τ˜ ) = −d(σ, χ) for τ˜ = (0, χ) ∈ X. Thus system (2.11) becomes: find σ˜ ∈ X and u˜ ∈ V such that ˜ = 0, τ˜ = (τ, 0) ∈ X, a1 (σ˜ , τ˜ ) + b1 (τ˜ , u) , (3.54) b(σ˜ , v˜ ) − c(ρu, v) = −f , v , ∀v˜ ∈ V ˜ = g, χ , ∀χ ∈ Q, c(κ p, χ) + b2 (τ˜ , u) where u, v = Ω uv dx. Similar for the discrete system (3.1). Indeed, letting τh = 0 in (3.1a) and recalling that Mh is the space consisting of piecewise constants, we get tr(σh )χh dx − χh λnh ds = 0, ∀χh ∈ Mh . (3.55) K∈Th∂K
Ω
Combining (3.55) and (3.1c), the discrete problem (3.1) can be reformulated as follows: µσ : τ + ν tr(σ ) tr(τ ) + div µτ + ν tr(τ )δ · uh h h h h h h Ω K K∈ T h − (τh )nt λt h ds = 0, µτh + ν tr(τh )δ nn λnh ds − µ ∂K ∂K K∈T K∈ T h h κ ph χh div U + χh λnh ds = gχh , − K Ω K∈ Th K∈Th ∂K div µσ + ν tr(σ )δ · v + ρ u · v div U = − f vh, h h h h h K Ω Ω K∈Th + ν tr(σ )δ η ds + ph ηnh ds = 0, µσ − h h nn nh ∂K ∂K K∈Th K∈Th (σh )nt ηt h ds = 0, −µ ∂K
(3.56)
K∈Th
for all (τ˜h = (τh , χh ); v˜ h = (v h , ηnh , ηt h )) ∈ X h × Mh . The algebraic equations generated by (3.56) are of the form σ 0 A 0 BT DT E T 0 HT 0 p G 0 C 0 0 u = F , B 0 Z D H 0 0 0 q 0 E 0 0 0 0 r 0
(3.57)
where σ, p, u, q and r denote by the degrees of freedom for σh , ph , uh , λnh and λt h , respectively. Since A is (symmetric) positive definite and block diagonal, the unknown σ can be eliminated, i.e., (3.58) σ = − A−1 B T u + A−1 D T q + A−1 E T r . Putting σ into the third and fourth equations of (3.57), we have
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
F = Z − BA−1 B T u − BA−1 D T q − BA−1 E T r, −1
−1
−1 T
DA B u + DA D q + DA E r − Hp = 0. T
T
291
(3.59) (3.60)
Set A∗ = Z − BA−1 B T . If | div U |∞ is small enough, then A∗ will be invertible and the unknown u is eliminated and −1 T −1 −1 T −1 u = A−1 ∗ BA D q + A∗ BA E r + A∗ F.
(3.61)
Putting (3.61) into (3.60), and inserting (3.61) into (3.58) and using the equation Eσ = 0 we have
q + Yr − Hp = F
1 , X T
2 , (3.62) Y q + Mr = F T H q − Cp = −G, where
−1 T T X = DA−1 B T A−1 , ∗ DA B −1 T
= X + DA D , X T −1 E + DA−1 E T , Y = D A−1 B T A−1 ∗ BA −1 T T + EA−1 E T , M = EA−1 B T A−1 ∗ EA B
1 = −DA−1 B T A−1 F ∗ F,
(3.63)
2 = −EA−1 B T A−1 F ∗ F.
in (3.63) is invertible, then (3.62) is written by If X
−1 H
−1 Y p Λ HT X −C + HT X 1 T −1 T = , T −1
r Λ M−Y X Y H X Y 2
(3.64)
2 − Y T X
−1 F
1 − g and Λ2 = F
−1 F
1 . where Λ1 = HT X
Conclusion The linearized steady-state barotropic compressible viscous Navier–Stokes equations (1.1) are a system of PDEs of mixed type; the momentum equations form an elliptic subsystem in velocity variable and the continuity equation is a hyperbolic equation in pressure. Hence it is a non-self adjoint boundary value problem. However, from the numerical approximation of the mixed finite element method (3.1) by the lowest order Raviart–Thomas finite elements we have seen that the resulted discrete equations in (3.56) generate a symmetric form (3.57), which can be transformed into (3.64) and also the approximate solution of the discrete problem (3.1) converges to the true solution optimally.
References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] D.N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: Implementation, postprocessing and error estimates, RAIRO Modél. Math. Anal. Numér. 19 (1985) 7–32.
292
J.R. Kweon / Applied Numerical Mathematics 45 (2003) 275–292
[3] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337–344. [4] D.N. Arnold, R.S. Falk, A new mixed formulation for elasticity, Numer. Math. 16 (1971) 322–333. [5] I. Babuška, W.C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (1978) 736–754. [6] H. Beirão da Veiga, An Lp -theory for the n-dimensional, stationary, compressible Navier–Stokes equations, and the incompressible limit for compressible fluids. The equilibrium solutions, Commun. Math. Phys. 109 (1987) 229–248. [7] H. Beirão Da Veiga, Existence results in Sobolev spaces for a stationary transport equation, Ricerche Mat. Suppl. 36 (1987) 173–184. [8] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [9] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1991. [10] M. Farhloul, Méthodes d’Éléments Finis Mixtes et des Volumes Finis, Thèse, Université Laval, Québec, 1991. [11] M. Farhloul, M. Fortin, A new mixed finite element for the Stokes and elasticity problems, SIAM J. Numer. Anal. 30 (1993) 971–990. [12] V. Girault, P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, in: Springer Ser. Comput. Math., Vol. 5, Springer, Berlin, 1986. [13] R.B. Kellogg, B. Liu, A finite element method for the compressible Stokes equations, SIAM J. Numer. Anal. 1 (1996) 780–788. [14] J.R. Kweon, A discontinuous Galerkin method for convection-dominated compressible viscous Navier–Stokes equations with inflow boundary condition, SIAM J. Numer. Anal. 38 (2000) 699–717. [15] J.R. Kweon, An optimal order convergence for a weak formulation of the compressible Stokes system with inflow boundary condition, Numer. Math. 86 (2000) 305–318. [16] J.R. Kweon, A mixed finite element method for a compressible Stokes problem with high Reynolds number, Appl. Numer. Math. 38 (2001) 87–103. [17] J.R. Kweon, R.B. Kellogg, Compressible Navier–Stokes equations in a bounded domain with inflow boundary bondition, SIAM J. Math. Anal. 28 (1997) 94–108. [18] P.D. Lax, R.S. Phillips, Local boundary conditions for dissipative symmetric linear differential operators, Commun. Pure Appl. Math. 18 (1960) 427–455. [19] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, in: Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere, New York, 1980. [20] P.A. Raviart, J.M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: Lecture Notes in Math., Vol. 606, Springer, Berlin, 1977, pp. 292–315. [21] M.A. Saad, Compressible Fluid Flow, 2nd Edition, Prentice-Hall, London, 1993.