A ROBUST NUMERICAL METHOD FOR STOKES EQUATIONS ...
A ROBUST NUMERICAL METHOD FOR STOKES EQUATIONS BASED ON DIVERGENCE-FREE H (DIV) FINITE ELEMENT METHODS JUNPING WANG∗ , YANQIU WANG† , AND XIU YE‡ Abstract. A computational method based on a divergence-free H(div) approach is presented for the Stokes equations in this article. This method is designed to find velocity approximation in an exact divergence-free subspace of the corresponding H(div) finite element space. That is, the continuity equation is strongly enforced a priori and the pressure is eliminated from the linear system in calculation. A strength of this approach is that the saddle-point problem for Stokes equations is reduced to a symmetric positive definite problem in a subspace for which basis functions are readily available. The resulting discrete system can then be solved by using existing sophisticated solvers. The aim of this article is to demonstrate the efficiency and robustness of H(div) finite element methods for Stokes equations. The results not only confirm the existing theoretical results, but also reveal additional advantages of the method in dealing with discontinuous boundary conditions. Key words. finite element methods, divergence-free, Stokes equations AMS subject classifications. Primary, 65N15, 65N30, 76D07; Secondary, 35B45, 35J50
1. Introduction. One of the challenges in solving Navier-Stokes equations is that the velocity and the pressure variables are coupled in a mixed system with a difficult saddle-point property. Recent study has resulted in several efficient methods (e.g., projection methods and Uzawa type iterative methods) to overcome this difficulty. In this article, we are concerned with a divergence-free approach which essentially decouples the variables by computing an approximate velocity solution of the Stokes equations in a divergence-free subspace, weakly or exactly. The main objective of this article is to demonstrate the efficiency and numerical robustness for a newly developed H(div) finite element methods [30, 31]. In standard finite element methods for Navier-Stokes equations, both the pressure and the velocity variables are approximated simultaneously [15] by using finite element functions satisfying a stability condition – known as the inf-sup condition. This method, known as primitive variable approach, results in a large saddle-point problem for which most existing numerical solvers are less effective and robust than for definite systems. While such saddle-point systems can be reduced to definite problems for the velocity unknown defined in weakly divergence-free subspaces, it is generally challenging, if not impossible, to construct computationally feasible bases for the resulting (weakly) divergence-free subspaces. This difficulty significantly limited the advantages of divergence-free approach in solving Stokes and Navier-Stokes equations. In addition to the primitive variable approach with standard Galerkin methods, attention was recently paid to numerical methods by using discontinuous finite ele∗ Division of Mathematical Sciences, National Science Foundation, Arlington, VA 22230 (jwang@ nsf.gov). This material was based on work supported by the National Science Foundation, while working at the Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. † Department of Mathematics, Oklahoma State University, Stillwater, OK 74075 ([email protected]). ‡ Department of Mathematics, University of Arkansas at Little Rock, Little Rock, AR 72204 ([email protected]). This research of Ye was supported in part by National Science Foundation Grant DMS-0612435
1
ments [6, 7, 16, 24, 25, 26, 34]. But, once again, most of them result in a numerical scheme in which the velocity is approximated in a weakly divergence-free subspace. Furthermore, the construction of computationally feasible basis functions was only possible for some special elements with certain particular orders [17, 18, 19, 21, 22, 32, 33, 35]. To the author’s knowledge, there is no systematic approach for constructing basis functions for (weakly) divergence-free finite element subspaces in the literature. An alternative way in approximating the Stokes equations is to use H(div) conforming finite elements [7, 30, 31]. The main motivation and advantage of using H(div) conforming elements for fluid flow problems is that the discrete velocity field will be globally exactly divergence-free, assuming that the fluid is incompressible. Since divergence-free exact functions can be written as the curl of a potential/stream function, divergence-free subspaces can then be constructed from taking curl of corresponding potential spaces. In two dimensional space, the potential functional space, also called stream functions, is well understood with basis functions readily available for computational purposes. In three dimensional spaces, vector potentials need to be considered and a large kernel of the curl operator adds to the complexity of the problem. In this paper, we focus on two-dimensional problems, with possible extension to three-dimensional cases. Another important feature of using H(div) conforming elements is that the method offers a more dynamic treatment of boundary conditions than the standard Garlekin methods. For example, in the H(div) method, the normal component of the velocity is set as an essential boundary condition and is strongly enforced, but the tangential component of the velocity is treated as a natural boundary condition and is weakly imposed. Thus, this approach empowers the H(div) method for problems with discontinuous boundary conditions. This nice feature is numerically illustrated in Section 5 for a lid driven cavity flow problem where the boundary condition is discontinuous. This paper is organized as follows. In section 2, we introduce some standard notation for Sobolev spaces. A variational formulation is also outlined for the Stokes problem. In section 3, we present a divergence-free H(div) finite element method by using the variational formulation developed previously. A detailed description of the divergence-free subspaces for H(div) conforming elements in two and three dimensions is also given. In Section 4, we present some numerical results for three test problems, each with a different configuration of boundary condition for the velocity. Finally in Section 5, we conduct some numerical tests for a lid driven cavity problem. 2. Preliminaries. Let us consider the Stokes equations with Dirichlet boundary conditions for the velocity variable. The problem seeks the velocity u and the pressure p in a certain functional spaces such that (2.1)
−∆u + ∇p = f
(2.2) (2.3)
∇·u=0 u=0
in Ω, in Ω, on ∂Ω,
where ∆, ∇, and ∇· denote the Laplacian, gradient, and divergence operators, respectively; Ω ⊂ Rd is the region occupied by the fluid; f = f (x) ∈ (L2 (Ω))d is the unit external volumetric force acting on the fluid at x ∈ Ω. A detailed discussion for inhomogeneous boundary value problems will be given in Section 3. For simplicity, the method will be presented for two-dimensional problems. Furthermore, we assume that Ω is a plane polygonal domain without cracks. 2
2.1. Notation. Let D be a bounded domain in R2 . We use standard definitions for the Sobolev spaces H s (D) and their associated inner-products (·, ·)s,D , norms k · ks,D , and semi-norms | · |s,D for s ≥ 0 [1]. The space H 0 (D) coincides with L2 (D), for which the norm and the inner product are denoted by k·kD and (·, ·)D , respectively. When D = Ω, we shall drop the subscript D in the norm and inner product notation. We also use L20 (Ω) to denote the subspace of L2 (Ω) consisting of functions with mean value zero. Throughout the paper, we adopt the convention that a bold-face character denotes a vector. Define H(div; Ω) as the space of vector-valued functions by H(div; Ω) = v : v ∈ (L2 (Ω))2 , ∇ · v ∈ L2 (Ω) ,
and with the norm
kvkH(div;Ω) = kvk2 + k∇ · vk2
12
.
Let K ⊂ Ω be a triangle or quadrilateral. For any smooth vector-valued functions w and v, it follows from the divergence theorem that Z Z ∂w · v ds, (2.4) (−∆w) · vdK = (∇w, ∇v)K − ∂n K K ∂K where ds represents the boundary element, nK is the outward normal direction on ∂K, and (∇w, ∇v)K =
2 Z X
i,j=1
K
∂wi ∂vi dK. ∂xj ∂xj
Let τK be the tangential direction to ∂K so that nK and τK form a right-hand coordinate system. It follows from the representation v = (v · nK )nK + (v · τK )τK that (2.5)
∂(w · nK ) ∂(w · τK ) ∂w ·v = (v · nK ) + (v · τK ). ∂nK ∂nK ∂nK
2.2. A variational formulation. To solve the Stokes system (2.1)-(2.3), A discontinuous Galerkin type formulation and finite element discretization was introduced in [30]. Following the idea of [30], and for the reader’s convenience, we outline the finite element formulations in this subsection. Let Th be a quasi-uniform triangulation of Ω with characteristic mesh size h. Define finite element spaces Vh and Wh for the velocity and pressure, respectively, by Vh = {v ∈ H(div; Ω) : v|K ∈ Vk (K)
∀K ∈ Th ; v · n|∂Ω = 0},
Wh = {q ∈ L20 (Ω) : q|K ∈ Wk (K) ∀K ∈ Th }, where n is the outward normal direction on the boundary of Ω, (Vk (K), Wk (K)) can be any existing H(div) conforming finite element pairs with order k ≥ 1. For example, the Raviart-Thomas element of order k (RTk ) [27] or the Brezzi-DouglasMarini element of order k (BDMk ) [3]. Three new H(div) conforming elements have 3
n2
K1
τ1 n1
τ2
K2
e
Fig. 2.1. Normal and tangential vectors for neighbouring triangles.
been obtained by Wang and Ye [30]. They were created specially to solve the Stokes and Navier-Stokes equations. In this paper, we focus on the RTk and BDMk elements. Both of them satisfy the discrete inf-sup condition (also known as the LBB condition) [4]. Details of these elements are skipped since they can be found in numerous sources. Multiplying equation (2.1) by any test function v ∈ Vh , then using integration by parts and equation (2.4), we get Z X ∂u (2.6) · v ds − (p, ∇ · v) = (f , v). (∇u, ∇v)K − ∂K ∂nK K∈Th
Since v ∈ Vh , its normal component v · nK is continuous across each interior edge. Therefore, it follows from (2.5) that X Z X Z ∂u ∂(u · τK ) · v ds = v · τK ds. ∂n ∂nK K ∂K ∂K K∈Th
By defining (∇h u, ∇h v) = into (2.6), we obtain (2.7)
K∈Th
P
K∈Th (∇u, ∇v)K
(∇h u, ∇h v) − (p, ∇ · v) −
X Z
∂K
K∈Th
and substituting the above equation ∂(u · τK ) v · τK ds = (f , v). ∂nK
We now reformulate the boundary integrals in (2.7). Let e be an interior edge shared by two elements K1 and K2 . Denote unit normal vectors n1 , n2 and tangential directions τ1 , τ2 , respectively, on e for K1 and K2 (as shown in Figure 2.1). Define the average {{ · }} and jump [[ · ]] on e for vector-valued functions w as follows: 1 (n1 · ∇(w · τ1 )|∂K1 + n2 · ∇(w · τ2 )|∂K2 ) , 2 [[w]] = w|∂K1 · τ1 + w|∂K2 · τ2 .
{{ε(w)}} =
For boundary edge e = ∂K1 ∩ ∂Ω, the above two operations must be modified by {{ε(w)}} = n1 · ∇(w · τ1 )|∂K1 ,
[[w]] = w|∂K1 · τ1 .
Let Eh be the set of all edges, including boundary edges, in Th . For sufficiently 3 smooth u (e.g., u ∈ H 2 +ǫ (Ω) for some ǫ > 0), it is not hard to see that X Z XZ ∂(u · τK ) {{ε(u)}}[[v]]ds. v · τK ds = ∂nK ∂K e K∈Th
e∈Eh
4
Substituting the above into (2.7), we get (2.8)
XZ
(∇h u, ∇h v) − (∇ · v, p) −
e∈Eh
{{ε(u)}}[[v]]ds = (f , v).
e
This gives the first equation in the variational form. For the second, testing (2.2) against any q ∈ Wh yields (2.9)
(∇ · u, q) = 0.
Finally, let V (h) = Vh + (H s (Ω) ∩ H01 (Ω))2 , with s > 23 . Denote by XZ (2.10) {{ε(u)}}[[v]]ds, ao (u, v) = (∇h u, ∇h v) − e∈Eh
(2.11)
e
b(v, q) = (∇ · v, q),
two bilinear forms on V (h) × V (h) and V (h) × L20 (Ω). With the conditions specified in this paper, it can be proved that the exact solution (u; p) of the Stokes problem in 2D belongs to V (h) for some s > 32 . Readers are referred to [8, 9, 20, 25] for details. As a result, it follows from (2.8) and (2.9) that the exact solution of the 2D Stokes problem satisfies the following variational equations: (2.12) (2.13)
ao (u, v) − b(v, p) = (f , v) b(u, q) = 0
∀v ∈ Vh , ∀q ∈ Wh .
3. A divergence-free finite element method. The goal of this section is to describe a divergence-free finite element scheme based on the weak formulation (2.12)–(2.13). The div-free scheme will be numerically investigated with various test problems. 3.1. Finite element discretization and error estimates. First, we introduce a symmetric and a skew symmetric bilinear forms on V (h) × V (h) as follows: XZ as (w, v) = ao (w, v) + αh−1 {ε(v)}}[[w]] ds, e [[w]][[v]] − { e∈Eh
ans (w, v) = ao (w, v) +
e
XZ
e∈Eh
αh−1 {ε(v)}}[[w]] ds. e [[w]][[v]] + {
e
where α > 0 is a parameter, and he is the length of the edge e. Then the discrete problem can be written as follows (see [30] for detail): Algorithm 1. Find (uh ; ph ) ∈ Vh × Wh such that (3.1) (3.2)
a(uh , v) − b(v, ph ) = (f , v) b(uh , q) = 0
∀v ∈ Vh , ∀q ∈ Wh ,
where a(v, w) = as (v, w) or a(v, w) = ans (v, w) Remark 1. R As is usual in the practice of discontinuous Galerkin methods, in a(·, ·), the term e αh−1 e [[w]][[v]] ds is designed to ensure stability of the formulation. 5
The parameter α should be large enough to guarantee a good convergence rate, and as small as possible in order to keep the condition number of the discrete system low. The well-posedness of Algorithm 1 comes from the well-known discrete inf-sup condition [4] and the following lemma [30]: Lemma 3.1. The symmetric bilinear form as (·, ·) is coercive for sufficiently large α, and the skew symmetric bilinear form ans (·, ·) is coercive for any α > 0. Error estimates for Algorithm 1 are also given in [30]. We first introduce two norms ||| · |||1 and ||| · ||| on V (h) as follows: 2
|||v|||1 = |v|21,h +
(3.3)
X
2 h−1 e k[[v]]ke ,
e∈Eh 2
(3.4)
|||v||| =
2 |||v|||1
+
X
he k{{ε(v)}}k2e ,
e∈Eh
where
|v|21,h
=
P
2 K∈Th |v|1,K
and
kvk2e
=
Z
v · vds.
e
Theorem 3.2. Let (u; p) be the solution of (2.1)–(2.3) and (uh ; ph ) ∈ Vh × Wh be obtained from (3.1)–(3.2). Assume α large enough if a(v, w) = as (v, w). Then, there exists a constant C independent of h such that (3.5)
|||u − uh ||| + kp − ph k ≤ Chk (kukk+1 + kpkk ).
Furthermore, if the Stokes problem has the (H 2 )2 × H 1 -regularity property, then (3.6)
ku − uh k ≤ Chk+1 (kukk+1 + kpkk )
provided that (u; p) ∈ (H k+1 (Ω))2 × H k (Ω) with k ≥ 1. 3.2. Divergence-free scheme. Algorithm 1, as described in the previous section, introduces a coupled saddle point system. In order to solve this system efficiently, next, we will focus on decoupling it by using the divergence-free finite element method. Define the divergence-free subspace Dh of Vh by Dh ≡ {v ∈ Vh ; b(v, q) = 0,
∀q ∈ Wh }.
By properties of the BDMk and RTk finite element spaces, we have ∇ · Vh = Wh [4]. Therefore, it is easy to see that Dh = {v ∈ Vh ; ∇ · v = 0}. In other words, Dh is globally exactly divergence-free. We point out that this is usually not true for H 1 conforming velocity discretizations, and functions in Dh will then only be weakly divergence-free. By taking the test function in Dh , the discrete formulation (3.1)-(3.2) can be reduced into the following divergence-free finite element scheme: Algorithm 2. Find uh ∈ Dh such that for all v ∈ Dh (3.7)
a(uh , v) = (f , v).
Problem (3.7) is symmetric positive definite if we choose a(v, w) = as (v, w), which brings great advantage in numerical simulation since it can be solved by the 6
efficient conjugate gradient method. Also, developing preconditioners for symmetric positive definite equations is considerably easier than for other systems. Next, we give a computable basis for the divergence-free subspace Dh by using the potential from Helmholtz decomposition. In two-dimension, a divergence-free vector v admits a potential function φ and −∂y φ . v = curlφ := ∂x φ The 2D potential φ is usually called the stream-function in literature. The Helmholtz decomposition also holds for most 2D H(div) conforming finite element spaces [4, 12, 13]. For the Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements, the following result is well-known [4, 12, 13]: Theorem 3.3. There exists a one-to-one map curl : Sh → Dh where the streamfunction space Sh is defined as following: 1. for triangular RT element of order k ≥ 0 or BDM element of order k ≥ 1, Sh = {φ ∈ H01 (Ω); φ|K ∈ Pk+1 (K), K ∈ Th }; 2. for rectangular RT element of order k ≥ 0, Sh = {φ ∈ H01 (Ω); φ|K ∈ Qk+1 (K), K ∈ Th }; where Pk (K) = span{xi y j | 0 ≤ i + j ≤ k} and Qk (K) = span{xi y j | 0 ≤ i, j ≤ k}. According to the above Theorem, one can simply take curl of the nodal basis of Pk+1 or Qk+1 conforming spaces, and derive a computable basis for Dh . In three-dimension, the divergence-free vector field can similarly be identified as the curl of a vector potential. Here curl is the usual vector curl which maps 3D vectors into 3D vectors. Let Vh be the Raviart-Thomas discretization of H(div) with order k and Sh be the Nedelec edge discretization of H(curl) with order k. The following result is well known [2]: Theorem 3.4. Let Dh be the globally divergence-free subspace of Vh , then Dh = curlSh . The three-dimensional case is significantly more complicated than the two-dimensional case since the 3D curl operator has a fairly large kernel containing all gradient vectors. Hence it is usually not practical to derive a basis of Dh from a basis of Sh in three dimension. Although there are some results in this direction [5, 28], they are either complicated or carrying many limitations. For three-dimensional problem, another possible approach is to solve the problem directly on curl{φh }, where {φh } is a basis for Sh . Notice that curl{φh } is linearly dependent and hence is not a basis for Dh . Discretization using a linearly dependent spanning set leads to a singular linear algebraic system. However, we know that many Krylov subspace iterative solvers can handle singular systems well as long as the righthand side and the initial guess are orthogonal to the null space of the matrix. With careful design, the problem may still be solvable using Krylov subspace solvers. 3.3. Inhomogeneous boundary conditions. In practical computation, many Stokes problems are imposed with inhomogeneous boundary conditions for the velocity. Here we generalize the divergence-free finite element method to the following 7
problem: −∆u + ∇p = f in Ω, ∇ · u = 0 in Ω, u=g
on ∂Ω.
The boundary data can be decomposed into two parts: u·n=g·n
and
u · τ = g · τ.
The normal component u · n will be imposed as an essential boundary condition. In other words, we seek for discrete solutions from the following finite element space V˜h = {v ∈ H(div; Ω) : v|K ∈ Vk (K)
∀K ∈ Th ; v · n|∂Ω = Ih (g · n)},
where Ih (g · n) is a suitable nodal value interpolation on ∂Ω, based on the degrees of freedom of the discrete space for v · n on ∂Ω. The tangential component u · τ will be treated as a natural boundary condition, that is, it will be imposed weakly through boundary integrals. Let EhB be the set of all boundary edges in Th , then Algorithm 1 should be modified as follows. Find (uh ; ph ) ∈ V˜h × Wh such that a(u, v) − b(v, p) = (f , v) +
X Z
B e∈Eh
e
(g · τ )(αh−1 e v·τ ∓
∂(v · τ ) ) ds ∂n
b(u, q) = 0
∀v ∈ Vh , ∀q ∈ Wh ,
where the minus is taken when using the symmetric form as (·, ·) and the plus for ans (·, ·). As usual in treating inhomogeneous boundary problems, the test space should still carry the homogeneous boundary condition, which means v is in Vh instead of V˜h . Similar changes need to be made when using the divergence-free scheme. The right-hand side in algorithm 2 should be modified consequently. Furthermore, since the computational basis of Dh is derived from Sh , we need to be careful when imposing the essential boundary condition u·n = g·n. Indeed, we impose the essential boundary condition on Sh . To this end, let u ∈ Dh and φ ∈ Sh satisfies u = curlφ. Then ∂φ =u·n=g·n ∂τ
on ∂Ω.
Therefore, one only needs to impose the following essential boundary condition on Sh and compute the basis accordingly: φ(x) = φ(x0 ) +
Z
x
g · n ds
x ∈ ∂Ω.
x0
Here x0 is an arbitrary point on ∂Ω and the integral is taken counter-clockwisely along ∂Ω. This boundary condition for φ is unique only up to a constant φ(x0 ). However, the constant will go away after taking curl. In other words, one is free to choose any x0 and φ(x0 ) in the implementation. 8
4. Numerical experiments for the Stokes equations. Numerical results for two dimensional Stokes equations are presented in this section. The divergence-free finite element scheme introduced in Algorithm 2 is used. The main objective here is to numerically examine the accuracy and efficiency of the H(div) scheme. For simplicity, a rectangular computational domain with uniform rectangular grids are used in this numerical study. The Raviart-Thomas finite element of order k = 1 (RT1 ) is employed in the finite element discretization. For the symmetric bilinear formulation as (·, ·), the discrete system is solved by using the conjugate gradient method. The discretization from the non-symmetric formulation ans (·, ·) is solved by using the GMRES method. In both cases, a relative residual of ε = 1.0e − 8 is used as the stopping criteria. Let u be the exact velocity and uh be its divergence-free finite element approximation obtained from Algorithm 2. The error is calculated by computing various norms or semi-norms of Ph u − uh , where Ph is the nodal value interpolation into Q2 conforming finite element spaces. This provides an accurate and effective method for computing the error under various norms. To be more precise, we introduce the following notations !1/2 X , (∇(Ph u − uh ), ∇(Ph u − uh ))K E1 = K∈Th
E2 =
X
h−1 e k[[(Ph u
−
e∈Eh
uh )]]k2e
!1/2
,
E3 =
X
he k{{ε(Ph u −
uh )}}k2e
e∈Eh
!1/2
.
Clearly, we have |Ph u − uh |21 = E12 ,
2
|||Ph u − uh |||1 = E12 + E22 ,
2
|||Ph u − uh ||| = E12 + E22 + E32 .
The domain is given by Ω = (0, 1) × (0, 1), which is partitioned into uniform rectangular grids along the x- and y- axes. We examine the error for both the symmetric formulation and the non-symmetric formulation with various values for the stabilization parameter α. In particular, we would like to see the sensitivity of each formulation against the parameter values of α. 4.1. Some numerical results for test problem 1. The test problem 1 is a Stokes equation with exact solution given as follows: −x(x − 1)(2y − 1) u = 2xy(x − 1)(y − 1) , p = 0. y(y − 1)(2x − 1) It is clear that homogeneous boundary condition is satisfied by this velocity. Table 4.1 Numerical performance for test problem 1, using the symmetric formulation as with parameter value α = 100.
mesh iteration 8×8 147 16 × 16 633 32 × 32 2214 64 × 64 8432 128 × 128 32252 asym. order O(hk ), k =
E1 1.41e-02 7.04e-03 3.52e-03 1.76e-03 8.82e-04 1.00
E2 2.03e-04 7.33e-05 2.63e-05 9.35e-06 3.32e-06 1.48 9
E3 1.95e-02 7.18e-03 2.60e-03 9.31e-04 3.31e-04 1.47
k · kL2 5.05e-04 1.27e-04 3.19e-05 7.99e-06 2.00e-06 1.99
k · kL∞ 1.16e-03 3.19e-04 8.43e-05 2.17e-05 5.52e-06 1.93
Table 4.2 Error information for test problem 1, using the symmetric formulation as with α = 10.
mesh iteration 8×8 43 16 × 16 195 32 × 32 894 64 × 64 3863 128 × 128 15977 asym. order O(hk ), k =
E1 1.41e-02 7.07e-03 3.53e-03 1.76e-03 8.83e-04 1.00
E2 2.19e-03 8.03e-04 2.89e-04 1.03e-04 3.66e-05 1.48
E3 2.13e-02 7.91e-03 2.86e-03 1.02e-03 3.65e-04 1.47
k · kL2 4.17e-04 1.14e-04 3.02e-05 7.77e-06 1.97e-06 1.93
k · kL∞ 8.93e-04 2.57e-04 7.58e-05 2.07e-05 5.42e-06 1.84
Fig. 4.1. Rate of convergence for test problem 1, using the symmetric formulation as with α = 10.
Error for test problem 1 E
−2
10
0
E
1
−3
error
10
E2 2
−4
L
−5
L O(h)
10
∞
10
2
O(h ) −2
−1
10
10
h
The numerical results in tables 4.1 and 4.2 show an O(h) convergence for the ||| · ||| norm and O(h2 ) convergence for the L2 norm, which agrees well with the theoretical results in Theorem 3.2. The rate of convergence is illustrated in figure 4.1. It also shows that E2 and E3 , which are related to the jump on internal edges, are usually of higher order accuracy than the discrete semi H 1 -norm E1 . This phenomena was not predicted by any existing convergence theory. Furthermore, the L∞ norm of the error seems to be of order O(h2 ) accuracy, though no theoretical proof can be seen in any existing literature. Table 4.3 Numerical performance for test problem 1, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/53 16 × 16 1/418 24 × 24 190/482 32 × 32 501/501 asym. order O(hk ), k =
E1 1.40e-02 7.04e-03 4.69e-03 3.52e-03 1.0004
E2 1.99e-04 7.18e-05 3.94e-05 2.57e-05 1.4765 10
E3 1.92e-02 7.06e-03 3.90e-03 2.55e-03 1.4548
k · kL2 4.77e-04 1.17e-04 5.21e-05 2.92e-05 2.0138
k · kL∞ 1.12e-03 3.08e-04 1.41e-04 8.12e-05 1.8936
Table 4.4 Numerical performance for test problem 1, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/44 16 × 16 1/205 24 × 24 2/19 32 × 32 116/52 asym. order O(hk ), k =
E1 1.41e-02 7.06e-03 4.70e-03 3.52e-03 0.9998
E2 1.82e-03 6.61e-04 3.63e-04 2.37e-04 1.4716
E3 1.85e-02 6.69e-03 3.66e-03 2.38e-03 1.4787
k · kL2 2.53e-04 6.01e-05 2.62e-05 1.46e-05 2.0591
k · kL∞ 8.15e-04 2.32e-04 1.07e-04 6.17e-05 1.8591
The numerical results in tables 4.3 and 4.4 were for the non-symmetric formulation with the bilinear form ans (·; ·). Like in the symmetric case, it shows an O(h) convergence in the |||·|||-norm (or energy norm) and O(h2 ) convergence in the L2 norm, which are all of optimal-order. It should be pointed out that no optimal-order error estimate was known for the velocity approximation in the L2 norm. The restarted GMRES solver was employed when solving the linear system resulted from the non-symmetric formulation. This solver restarts in every 500 steps. The “iteration” is expressed as a pair of numbers m/n where m is the restart rounds and n is the number of steps in the current restart round. The total iteration number should be 500(m − 1) + n. For example, 1/53 indicates that 53 iterations were performed, and 2/19 would give 519 iterations. The maximum iteration number is 501/501, which is 250, 501 iterations. The restarted GMRES uses much less computer memory than the original GMRES since the original GMRES needs to save all orthogonal vectors generated in the Lanczos process. However, the restarted GMRES may have a slower convergence rate than the original one. In addition, it is not clear to the authors how the convergence of the restarted GMRES is characterized by the condition number of the underlying matrix. 4.2. Some numerical results for test problem 2. The test problem 2 is a Stokes equation with exact solution given as follows:
u=
sin (2 π x) cos (2 π y) , − cos (2 π x) sin (2 π y)
p = x2 + y 2 .
It is not hard to see that the following natural boundary conditions are satisfied by this velocity: u · n|∂Ω = 0,
u · τ |∂Ω 6= 0.
Note that the second boundary condition (along the tangential direction) only indicates that a natural boundary condition should be imposed for this problem when using the H(div) finite element methods. One would need to compute u · τ |∂Ω from the representation of the velocity u in the numerical test. The error information and rate of convergence for the numerical scheme are illustrated in Tables 4.5 and 4.6 for the symmetric formulation, and Tables 4.7 and 4.8 for the non-symmetric formulation. 11
Table 4.5 Numerical performance for test problem 2, using the symmetric formulation as with α = 100.
mesh iteration 8×8 139 16 × 16 602 32 × 32 1780 64 × 64 5190 128 × 128 15064 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 1.40e-02 4.59e-03 1.59e-03 5.55e-04 1.96e-04 1.54
E3 1.21e+00 4.34e-01 1.54e-01 5.48e-02 1.94e-02 1.49
k · kL2 3.25e-02 7.97e-03 1.99e-03 4.98e-04 1.25e-04 2.00
k · kL∞ 6.41e-02 1.72e-02 4.49e-03 1.15e-03 2.92e-04 1.95
Table 4.6 Error information for test problem 2, using as , α = 10.
mesh iteration 8×8 32 16 × 16 159 32 × 32 751 64 × 64 3364 128 × 128 14305 asym. order O(hk ), k =
E1 1.01e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 1.40e-01 4.97e-02 1.74e-02 6.13e-03 2.15e-03 1.51
E3 1.33e+00 4.79e-01 1.70e-01 6.05e-02 2.14e-02 1.49
k · kL2 2.62e-02 7.11e-03 1.87e-03 4.83e-04 1.22e-04 1.94
k · kL∞ 5.25e-02 1.42e-02 4.07e-03 1.10e-03 2.87e-04 1.87
Table 4.7 Numerical performance for test problem 2, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/39 16 × 16 1/318 24 × 24 13/30 32 × 32 46/188 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.36e-01 2.51e-01 1.0074
E2 1.37e-02 4.49e-03 2.41e-03 1.55e-03 1.5763
E3 1.20e+00 4.27e-01 2.33e-01 1.52e-01 1.4908
k · kL2 3.16e-02 7.60e-03 3.35e-03 1.90e-03 2.0290
k · kL∞ 6.28e-02 1.67e-02 7.63e-03 4.35e-03 1.9244
Table 4.8 Numerical performance for test problem 2, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/33 16 × 16 1/167 24 × 24 1/427 32 × 32 9/36 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.36e-01 2.52e-01 1.0017
E2 1.18e-01 4.09e-02 2.21e-02 1.42e-02 1.5290
E3 1.21e+00 4.17e-01 2.24e-01 1.44e-01 1.5334
k · kL2 2.02e-02 4.62e-03 2.00e-03 1.11e-03 2.0938
k · kL∞ 4.04e-02 1.11e-02 5.35e-03 3.12e-03 1.8449
4.3. Some numerical results for test problem 3. In test problem 3, the exact solution of the Stokes equations is given by cos (2 π x) sin (2 π y) u= , p = x2 + y 2 . − sin (2 π x) cos (2 π y) 12
The following essential boundary conditions are satisfied by this velocity:
u · n|∂Ω 6= 0,
u · τ |∂Ω = 0.
Again, one would need to compute u · n|∂Ω from the representation of this velocity in doing numerical computation. The convergence and error profile for different mesh configuration are illustrated in Tables 4.9–4.12. Table 4.9 Numerical performance for test problem 3, using the symmetric formation as with α = 100.
mesh iter. 8×8 152 16 × 16 593 32 × 32 1921 64 × 64 8367 128 × 128 27917 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 4.16e-02 7.77e-03 1.42e-03 2.60e-04 4.79e-05 2.44
E3 4.29e-01 9.81e-02 2.35e-02 5.76e-03 1.43e-03 2.05
k · kL2 2.40e-02 5.73e-03 1.42e-03 3.52e-04 8.87e-05 2.02
k · kL∞ 5.24e-02 1.32e-02 3.38e-03 8.27e-04 2.10e-04 1.99
Table 4.10 Error information for test problem 3, using the symmetric formulation as with α = 10.
mesh iter. 8×8 48 16 × 16 218 32 × 32 951 64 × 64 3964 128 × 128 16122 asym. order O(hk ), k =
E1 1.00e+00 5.04e-01 2.51e-01 1.25e-01 6.31e-02 1.00
E2 7.44e-02 1.60e-02 3.61e-03 8.50e-04 2.05e-04 2.12
E3 4.66e-01 1.08e-01 2.60e-02 6.37e-03 1.58e-03 2.05
k · kL2 2.27e-02 5.65e-03 1.41e-03 3.52e-04 8.83e-05 2.00
k · kL∞ 5.08e-02 1.31e-02 3.29e-03 8.24e-04 2.08e-04 1.99
Table 4.11 Numerical performance for test problem 3, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/101 16 × 16 2/21 24 × 24 111/414 32 × 32 350/304 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 3.36e-01 2.52e-01 1.0095
E2 4.15e-02 7.76e-03 2.87e-03 1.42e-03 2.4340 13
E3 4.26e-01 9.76e-02 4.22e-02 2.35e-02 2.0928
k · kL2 2.34e-02 5.60e-03 2.48e-03 1.41e-03 2.0293
k · kL∞ 5.23e-02 1.29e-02 5.77e-03 3.30e-03 1.9965
Table 4.12 Numerical performance for test problem 3, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/50 16 × 16 1/236 24 × 24 3/265 32 × 32 40/290 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.38e-01 2.51e-01 1.0005
E2 6.38e-02 1.35e-02 2.60e-02 3.02e-03 1.7712
E3 4.37e-01 1.01e-01 5.04e-02 2.58e-02 2.0145
k · kL2 1.92e-02 4.74e-03 2.21e-03 1.51e-03 1.8642
k · kL∞ 4.90e-02 1.26e-02 5.67e-03 3.45e-03 1.9272
4.4. Condition number and error dependency on α. We also tested how the condition number of the discrete system and the error depend on various values of the stabilization parameter α. The existing theory predicted that, for the symmetric formulation, the numerical method is stable and accurate for sufficiently large values of α. Since the discrete system is symmetric and positive definite, the condition number of the discrete system can be conveniently calculated using estimates for extreme eigenvalues from the conjugate gradient (Lanczos-type) process. To this end, we solve the test problem 1 with different values of α, and compare the results in Tables 4.13 and 4.14. The condition number for the symmetric case seems to depend linearly on α. As α becomes larger, the tangential jump across internal edges will become smaller. Hence it is expected that E2 gets smaller as α increases. From Tables 4.13 and 4.14, it seems that the error in other norms are not affected much by α when α is large enough. It can also be seen that the non-symmetric formulation is not sensitive to the change of α, especially with small values. Table 4.13 Condition number and error for different values of α for test problem 1 on the 16 × 16 grid, with symmetric formulation as .
α 1 2 4 8 16 32 64 128
iteration 252 166 172 183 230 310 398 488
condition – 5.91e+04 8.78e+04 1.55e+05 2.97e+05 5.87e+05 1.17e+06 2.33e+06
E1 7.01e-02 1.19e-02 7.46e-03 7.10e-03 7.05e-03 7.04e-03 7.04e-03 7.04e-03
E2 6.74e-02 8.68e-03 2.44e-03 1.03e-03 4.82e-04 2.37e-04 1.15e-04 5.71e-05
E3 7.23e-02 1.76e-02 9.72e-03 8.15e-03 7.58e-03 7.34e-03 7.22e-03 7.17e-03
k · kL2 1.27e-03 2.43e-04 1.08e-04 1.12e-04 1.19e-04 1.24e-04 1.26e-04 1.27e-04
k · k∞ 3.28e-03 1.60e-03 4.74e-04 2.57e-04 2.80e-04 3.03e-04 3.15e-04 3.20e-04
5. Numerical experiments for a lid driven cavity problem. In this section, we report some numerical results on a lid driven cavity problem, for which the exact solution is not known. The 2D lid driven cavity problem describes the flow in a rectangular container which is driven by the uniform motion of one lid [29]. The lid driven cavity problem is one of the most popular benchmark problems for testing numerical schemes in fluid flow. One of the main difficulties of this problem is that it has a discontinuous velocity boundary condition and the standard primitive Galerkin methods have difficulty in dealing with such discontinuities without a further approximation of the boundary data. In two dimensional case, this boundary condition results in corner singularities for the solution. 14
Fig. 5.1. Streamline portrait of the lid driven cavity problem, obtained from the symmetric formulation as with various values of the stabilization parameter α. stream contour, alpha=1
stream contour, alpha=8
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
stream contour, alpha=128 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
stream contour, alpha=512
1
0 0
0.2
0 0
1
0.2
0.4
0.6
0.8
1
Fig. 5.2. Streamline portrait for the lid driven cavity problem, obtained from the non-symmetric, but absolutely stable formulation ans with various values of the stabilization parameter α. stream contour, alpha=0.01 (NS)
stream contour, alpha=0.1 (NS)
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
stream contour, alpha=1 (NS) 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
stream contour, alpha=8 (NS)
1
0 0
0.2
0 0
1
15
0.2
0.4
0.6
0.8
1
Table 4.14 Numerical results for different values of α for test problem 1 on the 16 × 16 grid, with nonsymmetric formulation ans .
α 1 2 4 8 16 32 64 128
iteration 1/187 1/182 1/183 1/195 1/237 1/303 1/374 1/440
E1 7.73e-03 7.35e-03 7.15e-03 7.07e-03 7.04e-03 7.04e-03 7.04e-03 7.04e-03
E2 4.18e-03 2.55e-03 1.47e-03 8.10e-04 4.27e-04 2.20e-04 1.11e-04 5.62e-05
E3 5.28e-03 5.76e-03 6.24e-03 6.60e-03 6.83e-03 6.96e-03 7.03e-03 7.07e-03
k · kL2 4.10e-04 2.55e-04 1.27e-04 6.18e-05 7.40e-05 9.76e-05 1.12e-04 1.20e-04
k · kL∞ 1.16e-03 7.64e-04 4.68e-04 2.75e-04 2.20e-04 2.71e-04 2.98e-04 3.12e-04
The 2D lid driven cavity problem was formulated in Ω = (0, 1)2 , with boundary condition u = (1, 0)t on the top lid and u = (0, 0)t elsewhere. In most direct numerical simulations of this problem, one has to choose explicitly the essential boundary data on two top corners. Popular choices are the “leaky” type, where u = (1, 0)t ; the “non-leaky” type, where u = (0, 0)t ; and the “regularized” type, where the boundary data on the top lid is replaced by a smooth function which vanishes at two corners. We emphasize that, in our divergence-free H(div) method, the boundary data discontinuity no longer cause any difficulty in the numerical scheme. This is so because, as pointed out earlier, only the normal component of the boundary data u · n (which equals zero and is thus continuous) is imposed as the essential boundary condition. The tangential component u · τ , which carries the discontinuity, is imposed weakly through boundary integrals. Therefore, the H(div) finite element method is a natural fit to the lid driven cavity problem (of cause the method was not motivated by the lid driven cavity problem). The computational results for the two-dimensional lid driven cavity problem are obtained by using a uniform 32 × 32 rectangular mesh with both symmetric and nonsymmetric finite element formulations. We are especially interested in seeing how the numerical solutions are affected by the change of α values. The streamline portraits are shown in Figures 5.1 and 5.2, and the velocity profiles are illustrated in Figures 5.3, 5.4, 5.5, 5.6, and 5.7. While we would like to leave readers to draw conclusions, we do like to point out the following obvious phenomena: • The symmetric scheme is non-stable when values of α are not sufficiently large. For example, the CG method did not converge for the linear system with α = 1. This means that the required positive definiteness of the linear system may fail to be valid for this case. • The non-symmetric scheme is stable regardless the value of α. Of course, the system’s coercivity gets weaker and weaker when the parameter α is getting closer to zero from positive. For small values of α, the continuity of the velocity approximation is less enforced. This can be seen from the numerical solution in Figures 5.5 and 5.6. But the discontinuity is suppressed when the value of α gets large as shown in Figure 5.7. • In the streamline portrait, the primary eddy and two corner eddies are clearly visible. • The non-symmetric finite element formulation is quite stable with respect to the parameter value of α. But extra care might be needed for solving 16
Fig. 5.3. The first velocity component profile for the lid driven cavity problem, obtained by using the symmetric formulation with various values of α. velocity u1, alpha=8
velocity u1, alpha=128
1
1
0.5
0.5
0
0
−0.5 1
−0.5 1
0 0.5
0 0.5
0.5 0
0.5 0
1
1
velocity u1, alpha=512
1 0.5 0 −0.5 1
0 0.5
0.5 0
1
the resulting non-symmetric matrix problem. GMRES was employed in our numerical experiments and the method has an acceptable convergence for the size of problem illustrated in this article. • The condition number for the resulting matrix is proportional to h−4 where h is the mesh size of the finite element partition. Therefore, neither CG nor GMRES provides a fast solver for the matrix problem. A good preconditioner is clearly needed in order to speedup the performance of CG and GMRES. Preconditioning issues will be further explored in forthcoming papers.
Fig. 5.6. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 0.1. velocity u1, alpha=0.1 (NS)
velocity u2, alpha=0.1 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5 0
0 0.5
0.5
0.5 0
1
17
1
Fig. 5.4. The second velocity component profile for the lid driven cavity problem, obtained by using the symmetric formulation with various values of α. velocity u2, alpha=8
velocity u2, alpha=128
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.4
−0.4
1
1
0 0.5
0 0.5
0.5 0
0.5 0
1
1
velocity u2, alpha=512
0.6 0.4 0.2 0 −0.2 −0.4 1
0 0.5
0.5 0
1
Fig. 5.5. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 0.01. velocity u1, alpha=0.01 (NS)
velocity u2, alpha=0.01 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5 0
0 0.5
0.5
0.5 0
1
1
Fig. 5.7. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 1. velocity u1, alpha=1 (NS)
velocity u2, alpha=1 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5
0.5 0
1
0 0.5
0.5
18 0
1
REFERENCES [1] R.A. Adams, J.J.F. Fournier, Sobolev Spaces, Elsevier, 2nd ed., Oxford, 2003. [2] Douglas N. Arnold, Richard S. Falk and Ragnar Winther, Multigrid in H(div) and H(curl), Numer. Math., 85 (2000) 197–218. [3] F. Brezzi, J. Douglas, and L. Marini, Two families of mixed finite elements for second order elliptic problem, Numer. Math., 47 (1985) 217–235. [4] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Elements, Springer-Verlag, New York, 1991. [5] Z. Cai, R. Parashkevov, T. Russell and X. Ye, Domain decomposition for a mixed finite element method in three dimensions, SIAM J. Numer. Anal., 41 (2003) 181–194. [6] B. Cockburn, G. Kanschat, and D. Schotzau, A locally conservative LDG method for the incompressible Navier-Stokes equations, Math. Comp., 74 (2005) 1067–1095. ¨ tzau, and C. Schwab, Local discontinuous Galerkin [7] B. Cockburn, G. Kanschat, D. Scho methods for the Stokes system, SIAM J. Numer. Anal., 40 (2002) 319–343. [8] M. Costabel and M. Dauge, Crack singularities for general elliptic systems, Math. Nachr., 235 (2002) 29–49. [9] M. Dauge, Elliptic Boundary Value Problems in Corner Domains - Smoothness and Asymptotics of Solutions, Lecture Notes in Math. 1341, Springer-Verlag, Berlin, 1988. [10] M. Dauge, Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. Part I. Linearized Equations, SIAM J. Math. Anal., 20 (1989) 74–97. [11] J. Douglas, Z. Cai and X. Ye, A Stable Quadrilateral nonconforming Element for the NavierStokes Equations, Calcolo, 36 (1999) 215–232. [12] R.E. Ewing and J. Wang, Analysis of the Schwarz algorithm for mixed finite element methods, R.A.I.R.O. Modelisation Mathematique Analyse Numerique, 26 (1992) 739–756. [13] R.E. Ewing and J. Wang, Analysis of multilevel decomposition iterative methods for mixed finite element methods, R.A.I.R.O. Mathematical Modeling and Numerical Analysis, 28 (1994) 377–398. [14] V. Girault and J.-L. Lions, Two-grid finite-element schemes for the steady Navier-Stokes problem in polyhedra, Port. Math. (N.S.), 58 (2001) 25–57. [15] V. Girault and P.A. Raviart, Finite Element Methods for the Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, 1986. `re, and M.F. Wheeler, A discontinuous Gelerkin method with noncon[16] V. Girault, B. Rivie forming domain decomposition for Stokes and Navier-Stokes problems, Math. Comp., 74 (2004) 53–84. [17] D. Griffiths, Finite element for incompressible flow, Math. Meth. in Appl. Sci., 1 (1979), 16–31. [18] D. Griffiths, The construction of approximately divergence-free finite element, The Mathematics of Finite Element an Its Applications III, Ed. J.R. Whiteman, Academic Press, 1979. [19] D. Griffiths, An approximately divergence-free 9-node velocity element for incompressible flows, Inter. J. Num. Meth. in Fluid, 1 (1981) 323–346. [20] P. Grisvard, Boundary Value Problems in Non-Smooth Domains, Pitman, London, 1985. [21] K. Gustafson and R. Hartman, Divergence-free basis for finite element schemes in hydrodynamics, SIAM J. Numer. Anal., 20 (1983) 697–721. [22] K. Gustafson and R. Hartman, Graph theory and fluid dynamics, SIAM J. Alg. Disc. Meth., 6 (1985) 643–656. [23] M. D. Gunzburger, Finite Element Methods for Viscous Incompressible Flows, A Guide to Theory, Practice and Algorithms, Academic, San Diego, 1989. [24] O. A. Karakashian and W. N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM J. Numerical Analysis, 35 (1998) 93–120. [25] J. R. Kweon, Corner singularity for compressible Navier-Stokes problem, in Proceedings of the International Congress of Mathematics, Higher Education Press, China, 2002. [26] J. Liu and C. Shu, A high order discontinuous Galerkin method for 2D incompressible flows, J. Comput. Phys., 160 (2000) 577–596. [27] P.A. Raviart and J.M. Thomas, A mixed finite element method for second order elliptic problems, Mathematical aspects of the finite element method, Eds. I. Galligani and E. Magenes, Lecture notes in Mathematics, Vol. 606., Springer-Verlag, 1977. [28] R. Scheichl, Decoupling three-dimensional mixed problems using divergence-free finite elements, SIAM J. Sci. Comp., 23 (2002) 1752–1776. [29] P.N. Shankar and M.D. Deshpande, Fluid Mechanics in the driven cavity, Annu. Rev. Fluid Mech., 32 (2000) 93–136. 19
[30] J. Wang and X. Ye, New finite element methods in computational fluid dynamics by H(div) elements, SIAM J. Numer. Anal., 45 (2007) 1269–1286. [31] J. Wang, X. Wang and X. Ye, Finite Element Methods for the Navier-Stokes equations by H(div) Elements, Journal of Computational Mathematics, 26 (2008) 410–436. [32] X. Ye and C. Hall, A Discrete divergence free basis for finite element methods, Numerical Algorithms, 16 (1997) 365–380. [33] X. Ye and C. Hall, The Construction of null basis for a discrete divergence operator, J. Computational and Applied Mathematics, 58 (1995) 117–133. [34] X. Ye, Discontinuous stable elements for the incompressible flow, Advances in Computational Mathematics, 20 (2004) 333–345. [35] X. Ye, Construction divergence free space for incompressible Navier-Stokes equations, Ph.D thesis, 1990.
20
1. Introduction. One of the challenges in solving Navier-Stokes equations is that the velocity and the pressure variables are coupled in a mixed system with a difficult saddle-point property. Recent study has resulted in several efficient methods (e.g., projection methods and Uzawa type iterative methods) to overcome this difficulty. In this article, we are concerned with a divergence-free approach which essentially decouples the variables by computing an approximate velocity solution of the Stokes equations in a divergence-free subspace, weakly or exactly. The main objective of this article is to demonstrate the efficiency and numerical robustness for a newly developed H(div) finite element methods [30, 31]. In standard finite element methods for Navier-Stokes equations, both the pressure and the velocity variables are approximated simultaneously [15] by using finite element functions satisfying a stability condition – known as the inf-sup condition. This method, known as primitive variable approach, results in a large saddle-point problem for which most existing numerical solvers are less effective and robust than for definite systems. While such saddle-point systems can be reduced to definite problems for the velocity unknown defined in weakly divergence-free subspaces, it is generally challenging, if not impossible, to construct computationally feasible bases for the resulting (weakly) divergence-free subspaces. This difficulty significantly limited the advantages of divergence-free approach in solving Stokes and Navier-Stokes equations. In addition to the primitive variable approach with standard Galerkin methods, attention was recently paid to numerical methods by using discontinuous finite ele∗ Division of Mathematical Sciences, National Science Foundation, Arlington, VA 22230 (jwang@ nsf.gov). This material was based on work supported by the National Science Foundation, while working at the Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation. † Department of Mathematics, Oklahoma State University, Stillwater, OK 74075 ([email protected]). ‡ Department of Mathematics, University of Arkansas at Little Rock, Little Rock, AR 72204 ([email protected]). This research of Ye was supported in part by National Science Foundation Grant DMS-0612435
1
ments [6, 7, 16, 24, 25, 26, 34]. But, once again, most of them result in a numerical scheme in which the velocity is approximated in a weakly divergence-free subspace. Furthermore, the construction of computationally feasible basis functions was only possible for some special elements with certain particular orders [17, 18, 19, 21, 22, 32, 33, 35]. To the author’s knowledge, there is no systematic approach for constructing basis functions for (weakly) divergence-free finite element subspaces in the literature. An alternative way in approximating the Stokes equations is to use H(div) conforming finite elements [7, 30, 31]. The main motivation and advantage of using H(div) conforming elements for fluid flow problems is that the discrete velocity field will be globally exactly divergence-free, assuming that the fluid is incompressible. Since divergence-free exact functions can be written as the curl of a potential/stream function, divergence-free subspaces can then be constructed from taking curl of corresponding potential spaces. In two dimensional space, the potential functional space, also called stream functions, is well understood with basis functions readily available for computational purposes. In three dimensional spaces, vector potentials need to be considered and a large kernel of the curl operator adds to the complexity of the problem. In this paper, we focus on two-dimensional problems, with possible extension to three-dimensional cases. Another important feature of using H(div) conforming elements is that the method offers a more dynamic treatment of boundary conditions than the standard Garlekin methods. For example, in the H(div) method, the normal component of the velocity is set as an essential boundary condition and is strongly enforced, but the tangential component of the velocity is treated as a natural boundary condition and is weakly imposed. Thus, this approach empowers the H(div) method for problems with discontinuous boundary conditions. This nice feature is numerically illustrated in Section 5 for a lid driven cavity flow problem where the boundary condition is discontinuous. This paper is organized as follows. In section 2, we introduce some standard notation for Sobolev spaces. A variational formulation is also outlined for the Stokes problem. In section 3, we present a divergence-free H(div) finite element method by using the variational formulation developed previously. A detailed description of the divergence-free subspaces for H(div) conforming elements in two and three dimensions is also given. In Section 4, we present some numerical results for three test problems, each with a different configuration of boundary condition for the velocity. Finally in Section 5, we conduct some numerical tests for a lid driven cavity problem. 2. Preliminaries. Let us consider the Stokes equations with Dirichlet boundary conditions for the velocity variable. The problem seeks the velocity u and the pressure p in a certain functional spaces such that (2.1)
−∆u + ∇p = f
(2.2) (2.3)
∇·u=0 u=0
in Ω, in Ω, on ∂Ω,
where ∆, ∇, and ∇· denote the Laplacian, gradient, and divergence operators, respectively; Ω ⊂ Rd is the region occupied by the fluid; f = f (x) ∈ (L2 (Ω))d is the unit external volumetric force acting on the fluid at x ∈ Ω. A detailed discussion for inhomogeneous boundary value problems will be given in Section 3. For simplicity, the method will be presented for two-dimensional problems. Furthermore, we assume that Ω is a plane polygonal domain without cracks. 2
2.1. Notation. Let D be a bounded domain in R2 . We use standard definitions for the Sobolev spaces H s (D) and their associated inner-products (·, ·)s,D , norms k · ks,D , and semi-norms | · |s,D for s ≥ 0 [1]. The space H 0 (D) coincides with L2 (D), for which the norm and the inner product are denoted by k·kD and (·, ·)D , respectively. When D = Ω, we shall drop the subscript D in the norm and inner product notation. We also use L20 (Ω) to denote the subspace of L2 (Ω) consisting of functions with mean value zero. Throughout the paper, we adopt the convention that a bold-face character denotes a vector. Define H(div; Ω) as the space of vector-valued functions by H(div; Ω) = v : v ∈ (L2 (Ω))2 , ∇ · v ∈ L2 (Ω) ,
and with the norm
kvkH(div;Ω) = kvk2 + k∇ · vk2
12
.
Let K ⊂ Ω be a triangle or quadrilateral. For any smooth vector-valued functions w and v, it follows from the divergence theorem that Z Z ∂w · v ds, (2.4) (−∆w) · vdK = (∇w, ∇v)K − ∂n K K ∂K where ds represents the boundary element, nK is the outward normal direction on ∂K, and (∇w, ∇v)K =
2 Z X
i,j=1
K
∂wi ∂vi dK. ∂xj ∂xj
Let τK be the tangential direction to ∂K so that nK and τK form a right-hand coordinate system. It follows from the representation v = (v · nK )nK + (v · τK )τK that (2.5)
∂(w · nK ) ∂(w · τK ) ∂w ·v = (v · nK ) + (v · τK ). ∂nK ∂nK ∂nK
2.2. A variational formulation. To solve the Stokes system (2.1)-(2.3), A discontinuous Galerkin type formulation and finite element discretization was introduced in [30]. Following the idea of [30], and for the reader’s convenience, we outline the finite element formulations in this subsection. Let Th be a quasi-uniform triangulation of Ω with characteristic mesh size h. Define finite element spaces Vh and Wh for the velocity and pressure, respectively, by Vh = {v ∈ H(div; Ω) : v|K ∈ Vk (K)
∀K ∈ Th ; v · n|∂Ω = 0},
Wh = {q ∈ L20 (Ω) : q|K ∈ Wk (K) ∀K ∈ Th }, where n is the outward normal direction on the boundary of Ω, (Vk (K), Wk (K)) can be any existing H(div) conforming finite element pairs with order k ≥ 1. For example, the Raviart-Thomas element of order k (RTk ) [27] or the Brezzi-DouglasMarini element of order k (BDMk ) [3]. Three new H(div) conforming elements have 3
n2
K1
τ1 n1
τ2
K2
e
Fig. 2.1. Normal and tangential vectors for neighbouring triangles.
been obtained by Wang and Ye [30]. They were created specially to solve the Stokes and Navier-Stokes equations. In this paper, we focus on the RTk and BDMk elements. Both of them satisfy the discrete inf-sup condition (also known as the LBB condition) [4]. Details of these elements are skipped since they can be found in numerous sources. Multiplying equation (2.1) by any test function v ∈ Vh , then using integration by parts and equation (2.4), we get Z X ∂u (2.6) · v ds − (p, ∇ · v) = (f , v). (∇u, ∇v)K − ∂K ∂nK K∈Th
Since v ∈ Vh , its normal component v · nK is continuous across each interior edge. Therefore, it follows from (2.5) that X Z X Z ∂u ∂(u · τK ) · v ds = v · τK ds. ∂n ∂nK K ∂K ∂K K∈Th
By defining (∇h u, ∇h v) = into (2.6), we obtain (2.7)
K∈Th
P
K∈Th (∇u, ∇v)K
(∇h u, ∇h v) − (p, ∇ · v) −
X Z
∂K
K∈Th
and substituting the above equation ∂(u · τK ) v · τK ds = (f , v). ∂nK
We now reformulate the boundary integrals in (2.7). Let e be an interior edge shared by two elements K1 and K2 . Denote unit normal vectors n1 , n2 and tangential directions τ1 , τ2 , respectively, on e for K1 and K2 (as shown in Figure 2.1). Define the average {{ · }} and jump [[ · ]] on e for vector-valued functions w as follows: 1 (n1 · ∇(w · τ1 )|∂K1 + n2 · ∇(w · τ2 )|∂K2 ) , 2 [[w]] = w|∂K1 · τ1 + w|∂K2 · τ2 .
{{ε(w)}} =
For boundary edge e = ∂K1 ∩ ∂Ω, the above two operations must be modified by {{ε(w)}} = n1 · ∇(w · τ1 )|∂K1 ,
[[w]] = w|∂K1 · τ1 .
Let Eh be the set of all edges, including boundary edges, in Th . For sufficiently 3 smooth u (e.g., u ∈ H 2 +ǫ (Ω) for some ǫ > 0), it is not hard to see that X Z XZ ∂(u · τK ) {{ε(u)}}[[v]]ds. v · τK ds = ∂nK ∂K e K∈Th
e∈Eh
4
Substituting the above into (2.7), we get (2.8)
XZ
(∇h u, ∇h v) − (∇ · v, p) −
e∈Eh
{{ε(u)}}[[v]]ds = (f , v).
e
This gives the first equation in the variational form. For the second, testing (2.2) against any q ∈ Wh yields (2.9)
(∇ · u, q) = 0.
Finally, let V (h) = Vh + (H s (Ω) ∩ H01 (Ω))2 , with s > 23 . Denote by XZ (2.10) {{ε(u)}}[[v]]ds, ao (u, v) = (∇h u, ∇h v) − e∈Eh
(2.11)
e
b(v, q) = (∇ · v, q),
two bilinear forms on V (h) × V (h) and V (h) × L20 (Ω). With the conditions specified in this paper, it can be proved that the exact solution (u; p) of the Stokes problem in 2D belongs to V (h) for some s > 32 . Readers are referred to [8, 9, 20, 25] for details. As a result, it follows from (2.8) and (2.9) that the exact solution of the 2D Stokes problem satisfies the following variational equations: (2.12) (2.13)
ao (u, v) − b(v, p) = (f , v) b(u, q) = 0
∀v ∈ Vh , ∀q ∈ Wh .
3. A divergence-free finite element method. The goal of this section is to describe a divergence-free finite element scheme based on the weak formulation (2.12)–(2.13). The div-free scheme will be numerically investigated with various test problems. 3.1. Finite element discretization and error estimates. First, we introduce a symmetric and a skew symmetric bilinear forms on V (h) × V (h) as follows: XZ as (w, v) = ao (w, v) + αh−1 {ε(v)}}[[w]] ds, e [[w]][[v]] − { e∈Eh
ans (w, v) = ao (w, v) +
e
XZ
e∈Eh
αh−1 {ε(v)}}[[w]] ds. e [[w]][[v]] + {
e
where α > 0 is a parameter, and he is the length of the edge e. Then the discrete problem can be written as follows (see [30] for detail): Algorithm 1. Find (uh ; ph ) ∈ Vh × Wh such that (3.1) (3.2)
a(uh , v) − b(v, ph ) = (f , v) b(uh , q) = 0
∀v ∈ Vh , ∀q ∈ Wh ,
where a(v, w) = as (v, w) or a(v, w) = ans (v, w) Remark 1. R As is usual in the practice of discontinuous Galerkin methods, in a(·, ·), the term e αh−1 e [[w]][[v]] ds is designed to ensure stability of the formulation. 5
The parameter α should be large enough to guarantee a good convergence rate, and as small as possible in order to keep the condition number of the discrete system low. The well-posedness of Algorithm 1 comes from the well-known discrete inf-sup condition [4] and the following lemma [30]: Lemma 3.1. The symmetric bilinear form as (·, ·) is coercive for sufficiently large α, and the skew symmetric bilinear form ans (·, ·) is coercive for any α > 0. Error estimates for Algorithm 1 are also given in [30]. We first introduce two norms ||| · |||1 and ||| · ||| on V (h) as follows: 2
|||v|||1 = |v|21,h +
(3.3)
X
2 h−1 e k[[v]]ke ,
e∈Eh 2
(3.4)
|||v||| =
2 |||v|||1
+
X
he k{{ε(v)}}k2e ,
e∈Eh
where
|v|21,h
=
P
2 K∈Th |v|1,K
and
kvk2e
=
Z
v · vds.
e
Theorem 3.2. Let (u; p) be the solution of (2.1)–(2.3) and (uh ; ph ) ∈ Vh × Wh be obtained from (3.1)–(3.2). Assume α large enough if a(v, w) = as (v, w). Then, there exists a constant C independent of h such that (3.5)
|||u − uh ||| + kp − ph k ≤ Chk (kukk+1 + kpkk ).
Furthermore, if the Stokes problem has the (H 2 )2 × H 1 -regularity property, then (3.6)
ku − uh k ≤ Chk+1 (kukk+1 + kpkk )
provided that (u; p) ∈ (H k+1 (Ω))2 × H k (Ω) with k ≥ 1. 3.2. Divergence-free scheme. Algorithm 1, as described in the previous section, introduces a coupled saddle point system. In order to solve this system efficiently, next, we will focus on decoupling it by using the divergence-free finite element method. Define the divergence-free subspace Dh of Vh by Dh ≡ {v ∈ Vh ; b(v, q) = 0,
∀q ∈ Wh }.
By properties of the BDMk and RTk finite element spaces, we have ∇ · Vh = Wh [4]. Therefore, it is easy to see that Dh = {v ∈ Vh ; ∇ · v = 0}. In other words, Dh is globally exactly divergence-free. We point out that this is usually not true for H 1 conforming velocity discretizations, and functions in Dh will then only be weakly divergence-free. By taking the test function in Dh , the discrete formulation (3.1)-(3.2) can be reduced into the following divergence-free finite element scheme: Algorithm 2. Find uh ∈ Dh such that for all v ∈ Dh (3.7)
a(uh , v) = (f , v).
Problem (3.7) is symmetric positive definite if we choose a(v, w) = as (v, w), which brings great advantage in numerical simulation since it can be solved by the 6
efficient conjugate gradient method. Also, developing preconditioners for symmetric positive definite equations is considerably easier than for other systems. Next, we give a computable basis for the divergence-free subspace Dh by using the potential from Helmholtz decomposition. In two-dimension, a divergence-free vector v admits a potential function φ and −∂y φ . v = curlφ := ∂x φ The 2D potential φ is usually called the stream-function in literature. The Helmholtz decomposition also holds for most 2D H(div) conforming finite element spaces [4, 12, 13]. For the Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements, the following result is well-known [4, 12, 13]: Theorem 3.3. There exists a one-to-one map curl : Sh → Dh where the streamfunction space Sh is defined as following: 1. for triangular RT element of order k ≥ 0 or BDM element of order k ≥ 1, Sh = {φ ∈ H01 (Ω); φ|K ∈ Pk+1 (K), K ∈ Th }; 2. for rectangular RT element of order k ≥ 0, Sh = {φ ∈ H01 (Ω); φ|K ∈ Qk+1 (K), K ∈ Th }; where Pk (K) = span{xi y j | 0 ≤ i + j ≤ k} and Qk (K) = span{xi y j | 0 ≤ i, j ≤ k}. According to the above Theorem, one can simply take curl of the nodal basis of Pk+1 or Qk+1 conforming spaces, and derive a computable basis for Dh . In three-dimension, the divergence-free vector field can similarly be identified as the curl of a vector potential. Here curl is the usual vector curl which maps 3D vectors into 3D vectors. Let Vh be the Raviart-Thomas discretization of H(div) with order k and Sh be the Nedelec edge discretization of H(curl) with order k. The following result is well known [2]: Theorem 3.4. Let Dh be the globally divergence-free subspace of Vh , then Dh = curlSh . The three-dimensional case is significantly more complicated than the two-dimensional case since the 3D curl operator has a fairly large kernel containing all gradient vectors. Hence it is usually not practical to derive a basis of Dh from a basis of Sh in three dimension. Although there are some results in this direction [5, 28], they are either complicated or carrying many limitations. For three-dimensional problem, another possible approach is to solve the problem directly on curl{φh }, where {φh } is a basis for Sh . Notice that curl{φh } is linearly dependent and hence is not a basis for Dh . Discretization using a linearly dependent spanning set leads to a singular linear algebraic system. However, we know that many Krylov subspace iterative solvers can handle singular systems well as long as the righthand side and the initial guess are orthogonal to the null space of the matrix. With careful design, the problem may still be solvable using Krylov subspace solvers. 3.3. Inhomogeneous boundary conditions. In practical computation, many Stokes problems are imposed with inhomogeneous boundary conditions for the velocity. Here we generalize the divergence-free finite element method to the following 7
problem: −∆u + ∇p = f in Ω, ∇ · u = 0 in Ω, u=g
on ∂Ω.
The boundary data can be decomposed into two parts: u·n=g·n
and
u · τ = g · τ.
The normal component u · n will be imposed as an essential boundary condition. In other words, we seek for discrete solutions from the following finite element space V˜h = {v ∈ H(div; Ω) : v|K ∈ Vk (K)
∀K ∈ Th ; v · n|∂Ω = Ih (g · n)},
where Ih (g · n) is a suitable nodal value interpolation on ∂Ω, based on the degrees of freedom of the discrete space for v · n on ∂Ω. The tangential component u · τ will be treated as a natural boundary condition, that is, it will be imposed weakly through boundary integrals. Let EhB be the set of all boundary edges in Th , then Algorithm 1 should be modified as follows. Find (uh ; ph ) ∈ V˜h × Wh such that a(u, v) − b(v, p) = (f , v) +
X Z
B e∈Eh
e
(g · τ )(αh−1 e v·τ ∓
∂(v · τ ) ) ds ∂n
b(u, q) = 0
∀v ∈ Vh , ∀q ∈ Wh ,
where the minus is taken when using the symmetric form as (·, ·) and the plus for ans (·, ·). As usual in treating inhomogeneous boundary problems, the test space should still carry the homogeneous boundary condition, which means v is in Vh instead of V˜h . Similar changes need to be made when using the divergence-free scheme. The right-hand side in algorithm 2 should be modified consequently. Furthermore, since the computational basis of Dh is derived from Sh , we need to be careful when imposing the essential boundary condition u·n = g·n. Indeed, we impose the essential boundary condition on Sh . To this end, let u ∈ Dh and φ ∈ Sh satisfies u = curlφ. Then ∂φ =u·n=g·n ∂τ
on ∂Ω.
Therefore, one only needs to impose the following essential boundary condition on Sh and compute the basis accordingly: φ(x) = φ(x0 ) +
Z
x
g · n ds
x ∈ ∂Ω.
x0
Here x0 is an arbitrary point on ∂Ω and the integral is taken counter-clockwisely along ∂Ω. This boundary condition for φ is unique only up to a constant φ(x0 ). However, the constant will go away after taking curl. In other words, one is free to choose any x0 and φ(x0 ) in the implementation. 8
4. Numerical experiments for the Stokes equations. Numerical results for two dimensional Stokes equations are presented in this section. The divergence-free finite element scheme introduced in Algorithm 2 is used. The main objective here is to numerically examine the accuracy and efficiency of the H(div) scheme. For simplicity, a rectangular computational domain with uniform rectangular grids are used in this numerical study. The Raviart-Thomas finite element of order k = 1 (RT1 ) is employed in the finite element discretization. For the symmetric bilinear formulation as (·, ·), the discrete system is solved by using the conjugate gradient method. The discretization from the non-symmetric formulation ans (·, ·) is solved by using the GMRES method. In both cases, a relative residual of ε = 1.0e − 8 is used as the stopping criteria. Let u be the exact velocity and uh be its divergence-free finite element approximation obtained from Algorithm 2. The error is calculated by computing various norms or semi-norms of Ph u − uh , where Ph is the nodal value interpolation into Q2 conforming finite element spaces. This provides an accurate and effective method for computing the error under various norms. To be more precise, we introduce the following notations !1/2 X , (∇(Ph u − uh ), ∇(Ph u − uh ))K E1 = K∈Th
E2 =
X
h−1 e k[[(Ph u
−
e∈Eh
uh )]]k2e
!1/2
,
E3 =
X
he k{{ε(Ph u −
uh )}}k2e
e∈Eh
!1/2
.
Clearly, we have |Ph u − uh |21 = E12 ,
2
|||Ph u − uh |||1 = E12 + E22 ,
2
|||Ph u − uh ||| = E12 + E22 + E32 .
The domain is given by Ω = (0, 1) × (0, 1), which is partitioned into uniform rectangular grids along the x- and y- axes. We examine the error for both the symmetric formulation and the non-symmetric formulation with various values for the stabilization parameter α. In particular, we would like to see the sensitivity of each formulation against the parameter values of α. 4.1. Some numerical results for test problem 1. The test problem 1 is a Stokes equation with exact solution given as follows: −x(x − 1)(2y − 1) u = 2xy(x − 1)(y − 1) , p = 0. y(y − 1)(2x − 1) It is clear that homogeneous boundary condition is satisfied by this velocity. Table 4.1 Numerical performance for test problem 1, using the symmetric formulation as with parameter value α = 100.
mesh iteration 8×8 147 16 × 16 633 32 × 32 2214 64 × 64 8432 128 × 128 32252 asym. order O(hk ), k =
E1 1.41e-02 7.04e-03 3.52e-03 1.76e-03 8.82e-04 1.00
E2 2.03e-04 7.33e-05 2.63e-05 9.35e-06 3.32e-06 1.48 9
E3 1.95e-02 7.18e-03 2.60e-03 9.31e-04 3.31e-04 1.47
k · kL2 5.05e-04 1.27e-04 3.19e-05 7.99e-06 2.00e-06 1.99
k · kL∞ 1.16e-03 3.19e-04 8.43e-05 2.17e-05 5.52e-06 1.93
Table 4.2 Error information for test problem 1, using the symmetric formulation as with α = 10.
mesh iteration 8×8 43 16 × 16 195 32 × 32 894 64 × 64 3863 128 × 128 15977 asym. order O(hk ), k =
E1 1.41e-02 7.07e-03 3.53e-03 1.76e-03 8.83e-04 1.00
E2 2.19e-03 8.03e-04 2.89e-04 1.03e-04 3.66e-05 1.48
E3 2.13e-02 7.91e-03 2.86e-03 1.02e-03 3.65e-04 1.47
k · kL2 4.17e-04 1.14e-04 3.02e-05 7.77e-06 1.97e-06 1.93
k · kL∞ 8.93e-04 2.57e-04 7.58e-05 2.07e-05 5.42e-06 1.84
Fig. 4.1. Rate of convergence for test problem 1, using the symmetric formulation as with α = 10.
Error for test problem 1 E
−2
10
0
E
1
−3
error
10
E2 2
−4
L
−5
L O(h)
10
∞
10
2
O(h ) −2
−1
10
10
h
The numerical results in tables 4.1 and 4.2 show an O(h) convergence for the ||| · ||| norm and O(h2 ) convergence for the L2 norm, which agrees well with the theoretical results in Theorem 3.2. The rate of convergence is illustrated in figure 4.1. It also shows that E2 and E3 , which are related to the jump on internal edges, are usually of higher order accuracy than the discrete semi H 1 -norm E1 . This phenomena was not predicted by any existing convergence theory. Furthermore, the L∞ norm of the error seems to be of order O(h2 ) accuracy, though no theoretical proof can be seen in any existing literature. Table 4.3 Numerical performance for test problem 1, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/53 16 × 16 1/418 24 × 24 190/482 32 × 32 501/501 asym. order O(hk ), k =
E1 1.40e-02 7.04e-03 4.69e-03 3.52e-03 1.0004
E2 1.99e-04 7.18e-05 3.94e-05 2.57e-05 1.4765 10
E3 1.92e-02 7.06e-03 3.90e-03 2.55e-03 1.4548
k · kL2 4.77e-04 1.17e-04 5.21e-05 2.92e-05 2.0138
k · kL∞ 1.12e-03 3.08e-04 1.41e-04 8.12e-05 1.8936
Table 4.4 Numerical performance for test problem 1, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/44 16 × 16 1/205 24 × 24 2/19 32 × 32 116/52 asym. order O(hk ), k =
E1 1.41e-02 7.06e-03 4.70e-03 3.52e-03 0.9998
E2 1.82e-03 6.61e-04 3.63e-04 2.37e-04 1.4716
E3 1.85e-02 6.69e-03 3.66e-03 2.38e-03 1.4787
k · kL2 2.53e-04 6.01e-05 2.62e-05 1.46e-05 2.0591
k · kL∞ 8.15e-04 2.32e-04 1.07e-04 6.17e-05 1.8591
The numerical results in tables 4.3 and 4.4 were for the non-symmetric formulation with the bilinear form ans (·; ·). Like in the symmetric case, it shows an O(h) convergence in the |||·|||-norm (or energy norm) and O(h2 ) convergence in the L2 norm, which are all of optimal-order. It should be pointed out that no optimal-order error estimate was known for the velocity approximation in the L2 norm. The restarted GMRES solver was employed when solving the linear system resulted from the non-symmetric formulation. This solver restarts in every 500 steps. The “iteration” is expressed as a pair of numbers m/n where m is the restart rounds and n is the number of steps in the current restart round. The total iteration number should be 500(m − 1) + n. For example, 1/53 indicates that 53 iterations were performed, and 2/19 would give 519 iterations. The maximum iteration number is 501/501, which is 250, 501 iterations. The restarted GMRES uses much less computer memory than the original GMRES since the original GMRES needs to save all orthogonal vectors generated in the Lanczos process. However, the restarted GMRES may have a slower convergence rate than the original one. In addition, it is not clear to the authors how the convergence of the restarted GMRES is characterized by the condition number of the underlying matrix. 4.2. Some numerical results for test problem 2. The test problem 2 is a Stokes equation with exact solution given as follows:
u=
sin (2 π x) cos (2 π y) , − cos (2 π x) sin (2 π y)
p = x2 + y 2 .
It is not hard to see that the following natural boundary conditions are satisfied by this velocity: u · n|∂Ω = 0,
u · τ |∂Ω 6= 0.
Note that the second boundary condition (along the tangential direction) only indicates that a natural boundary condition should be imposed for this problem when using the H(div) finite element methods. One would need to compute u · τ |∂Ω from the representation of the velocity u in the numerical test. The error information and rate of convergence for the numerical scheme are illustrated in Tables 4.5 and 4.6 for the symmetric formulation, and Tables 4.7 and 4.8 for the non-symmetric formulation. 11
Table 4.5 Numerical performance for test problem 2, using the symmetric formulation as with α = 100.
mesh iteration 8×8 139 16 × 16 602 32 × 32 1780 64 × 64 5190 128 × 128 15064 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 1.40e-02 4.59e-03 1.59e-03 5.55e-04 1.96e-04 1.54
E3 1.21e+00 4.34e-01 1.54e-01 5.48e-02 1.94e-02 1.49
k · kL2 3.25e-02 7.97e-03 1.99e-03 4.98e-04 1.25e-04 2.00
k · kL∞ 6.41e-02 1.72e-02 4.49e-03 1.15e-03 2.92e-04 1.95
Table 4.6 Error information for test problem 2, using as , α = 10.
mesh iteration 8×8 32 16 × 16 159 32 × 32 751 64 × 64 3364 128 × 128 14305 asym. order O(hk ), k =
E1 1.01e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 1.40e-01 4.97e-02 1.74e-02 6.13e-03 2.15e-03 1.51
E3 1.33e+00 4.79e-01 1.70e-01 6.05e-02 2.14e-02 1.49
k · kL2 2.62e-02 7.11e-03 1.87e-03 4.83e-04 1.22e-04 1.94
k · kL∞ 5.25e-02 1.42e-02 4.07e-03 1.10e-03 2.87e-04 1.87
Table 4.7 Numerical performance for test problem 2, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/39 16 × 16 1/318 24 × 24 13/30 32 × 32 46/188 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.36e-01 2.51e-01 1.0074
E2 1.37e-02 4.49e-03 2.41e-03 1.55e-03 1.5763
E3 1.20e+00 4.27e-01 2.33e-01 1.52e-01 1.4908
k · kL2 3.16e-02 7.60e-03 3.35e-03 1.90e-03 2.0290
k · kL∞ 6.28e-02 1.67e-02 7.63e-03 4.35e-03 1.9244
Table 4.8 Numerical performance for test problem 2, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/33 16 × 16 1/167 24 × 24 1/427 32 × 32 9/36 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.36e-01 2.52e-01 1.0017
E2 1.18e-01 4.09e-02 2.21e-02 1.42e-02 1.5290
E3 1.21e+00 4.17e-01 2.24e-01 1.44e-01 1.5334
k · kL2 2.02e-02 4.62e-03 2.00e-03 1.11e-03 2.0938
k · kL∞ 4.04e-02 1.11e-02 5.35e-03 3.12e-03 1.8449
4.3. Some numerical results for test problem 3. In test problem 3, the exact solution of the Stokes equations is given by cos (2 π x) sin (2 π y) u= , p = x2 + y 2 . − sin (2 π x) cos (2 π y) 12
The following essential boundary conditions are satisfied by this velocity:
u · n|∂Ω 6= 0,
u · τ |∂Ω = 0.
Again, one would need to compute u · n|∂Ω from the representation of this velocity in doing numerical computation. The convergence and error profile for different mesh configuration are illustrated in Tables 4.9–4.12. Table 4.9 Numerical performance for test problem 3, using the symmetric formation as with α = 100.
mesh iter. 8×8 152 16 × 16 593 32 × 32 1921 64 × 64 8367 128 × 128 27917 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 2.52e-01 1.26e-01 6.31e-02 1.00
E2 4.16e-02 7.77e-03 1.42e-03 2.60e-04 4.79e-05 2.44
E3 4.29e-01 9.81e-02 2.35e-02 5.76e-03 1.43e-03 2.05
k · kL2 2.40e-02 5.73e-03 1.42e-03 3.52e-04 8.87e-05 2.02
k · kL∞ 5.24e-02 1.32e-02 3.38e-03 8.27e-04 2.10e-04 1.99
Table 4.10 Error information for test problem 3, using the symmetric formulation as with α = 10.
mesh iter. 8×8 48 16 × 16 218 32 × 32 951 64 × 64 3964 128 × 128 16122 asym. order O(hk ), k =
E1 1.00e+00 5.04e-01 2.51e-01 1.25e-01 6.31e-02 1.00
E2 7.44e-02 1.60e-02 3.61e-03 8.50e-04 2.05e-04 2.12
E3 4.66e-01 1.08e-01 2.60e-02 6.37e-03 1.58e-03 2.05
k · kL2 2.27e-02 5.65e-03 1.41e-03 3.52e-04 8.83e-05 2.00
k · kL∞ 5.08e-02 1.31e-02 3.29e-03 8.24e-04 2.08e-04 1.99
Table 4.11 Numerical performance for test problem 3, using the non-symmetric formulation ans , with parameter value α = 100.
mesh iteration 8×8 1/101 16 × 16 2/21 24 × 24 111/414 32 × 32 350/304 asym. order O(hk ), k =
E1 1.02e+00 5.05e-01 3.36e-01 2.52e-01 1.0095
E2 4.15e-02 7.76e-03 2.87e-03 1.42e-03 2.4340 13
E3 4.26e-01 9.76e-02 4.22e-02 2.35e-02 2.0928
k · kL2 2.34e-02 5.60e-03 2.48e-03 1.41e-03 2.0293
k · kL∞ 5.23e-02 1.29e-02 5.77e-03 3.30e-03 1.9965
Table 4.12 Numerical performance for test problem 3, using the non-symmetric formulation ans , with parameter value α = 10.
mesh iteration 8×8 1/50 16 × 16 1/236 24 × 24 3/265 32 × 32 40/290 asym. order O(hk ), k =
E1 1.01e+00 5.04e-01 3.38e-01 2.51e-01 1.0005
E2 6.38e-02 1.35e-02 2.60e-02 3.02e-03 1.7712
E3 4.37e-01 1.01e-01 5.04e-02 2.58e-02 2.0145
k · kL2 1.92e-02 4.74e-03 2.21e-03 1.51e-03 1.8642
k · kL∞ 4.90e-02 1.26e-02 5.67e-03 3.45e-03 1.9272
4.4. Condition number and error dependency on α. We also tested how the condition number of the discrete system and the error depend on various values of the stabilization parameter α. The existing theory predicted that, for the symmetric formulation, the numerical method is stable and accurate for sufficiently large values of α. Since the discrete system is symmetric and positive definite, the condition number of the discrete system can be conveniently calculated using estimates for extreme eigenvalues from the conjugate gradient (Lanczos-type) process. To this end, we solve the test problem 1 with different values of α, and compare the results in Tables 4.13 and 4.14. The condition number for the symmetric case seems to depend linearly on α. As α becomes larger, the tangential jump across internal edges will become smaller. Hence it is expected that E2 gets smaller as α increases. From Tables 4.13 and 4.14, it seems that the error in other norms are not affected much by α when α is large enough. It can also be seen that the non-symmetric formulation is not sensitive to the change of α, especially with small values. Table 4.13 Condition number and error for different values of α for test problem 1 on the 16 × 16 grid, with symmetric formulation as .
α 1 2 4 8 16 32 64 128
iteration 252 166 172 183 230 310 398 488
condition – 5.91e+04 8.78e+04 1.55e+05 2.97e+05 5.87e+05 1.17e+06 2.33e+06
E1 7.01e-02 1.19e-02 7.46e-03 7.10e-03 7.05e-03 7.04e-03 7.04e-03 7.04e-03
E2 6.74e-02 8.68e-03 2.44e-03 1.03e-03 4.82e-04 2.37e-04 1.15e-04 5.71e-05
E3 7.23e-02 1.76e-02 9.72e-03 8.15e-03 7.58e-03 7.34e-03 7.22e-03 7.17e-03
k · kL2 1.27e-03 2.43e-04 1.08e-04 1.12e-04 1.19e-04 1.24e-04 1.26e-04 1.27e-04
k · k∞ 3.28e-03 1.60e-03 4.74e-04 2.57e-04 2.80e-04 3.03e-04 3.15e-04 3.20e-04
5. Numerical experiments for a lid driven cavity problem. In this section, we report some numerical results on a lid driven cavity problem, for which the exact solution is not known. The 2D lid driven cavity problem describes the flow in a rectangular container which is driven by the uniform motion of one lid [29]. The lid driven cavity problem is one of the most popular benchmark problems for testing numerical schemes in fluid flow. One of the main difficulties of this problem is that it has a discontinuous velocity boundary condition and the standard primitive Galerkin methods have difficulty in dealing with such discontinuities without a further approximation of the boundary data. In two dimensional case, this boundary condition results in corner singularities for the solution. 14
Fig. 5.1. Streamline portrait of the lid driven cavity problem, obtained from the symmetric formulation as with various values of the stabilization parameter α. stream contour, alpha=1
stream contour, alpha=8
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
stream contour, alpha=128 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
stream contour, alpha=512
1
0 0
0.2
0 0
1
0.2
0.4
0.6
0.8
1
Fig. 5.2. Streamline portrait for the lid driven cavity problem, obtained from the non-symmetric, but absolutely stable formulation ans with various values of the stabilization parameter α. stream contour, alpha=0.01 (NS)
stream contour, alpha=0.1 (NS)
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0 0
0.2
0.4
0.6
0.8
0 0
1
stream contour, alpha=1 (NS) 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
stream contour, alpha=8 (NS)
1
0 0
0.2
0 0
1
15
0.2
0.4
0.6
0.8
1
Table 4.14 Numerical results for different values of α for test problem 1 on the 16 × 16 grid, with nonsymmetric formulation ans .
α 1 2 4 8 16 32 64 128
iteration 1/187 1/182 1/183 1/195 1/237 1/303 1/374 1/440
E1 7.73e-03 7.35e-03 7.15e-03 7.07e-03 7.04e-03 7.04e-03 7.04e-03 7.04e-03
E2 4.18e-03 2.55e-03 1.47e-03 8.10e-04 4.27e-04 2.20e-04 1.11e-04 5.62e-05
E3 5.28e-03 5.76e-03 6.24e-03 6.60e-03 6.83e-03 6.96e-03 7.03e-03 7.07e-03
k · kL2 4.10e-04 2.55e-04 1.27e-04 6.18e-05 7.40e-05 9.76e-05 1.12e-04 1.20e-04
k · kL∞ 1.16e-03 7.64e-04 4.68e-04 2.75e-04 2.20e-04 2.71e-04 2.98e-04 3.12e-04
The 2D lid driven cavity problem was formulated in Ω = (0, 1)2 , with boundary condition u = (1, 0)t on the top lid and u = (0, 0)t elsewhere. In most direct numerical simulations of this problem, one has to choose explicitly the essential boundary data on two top corners. Popular choices are the “leaky” type, where u = (1, 0)t ; the “non-leaky” type, where u = (0, 0)t ; and the “regularized” type, where the boundary data on the top lid is replaced by a smooth function which vanishes at two corners. We emphasize that, in our divergence-free H(div) method, the boundary data discontinuity no longer cause any difficulty in the numerical scheme. This is so because, as pointed out earlier, only the normal component of the boundary data u · n (which equals zero and is thus continuous) is imposed as the essential boundary condition. The tangential component u · τ , which carries the discontinuity, is imposed weakly through boundary integrals. Therefore, the H(div) finite element method is a natural fit to the lid driven cavity problem (of cause the method was not motivated by the lid driven cavity problem). The computational results for the two-dimensional lid driven cavity problem are obtained by using a uniform 32 × 32 rectangular mesh with both symmetric and nonsymmetric finite element formulations. We are especially interested in seeing how the numerical solutions are affected by the change of α values. The streamline portraits are shown in Figures 5.1 and 5.2, and the velocity profiles are illustrated in Figures 5.3, 5.4, 5.5, 5.6, and 5.7. While we would like to leave readers to draw conclusions, we do like to point out the following obvious phenomena: • The symmetric scheme is non-stable when values of α are not sufficiently large. For example, the CG method did not converge for the linear system with α = 1. This means that the required positive definiteness of the linear system may fail to be valid for this case. • The non-symmetric scheme is stable regardless the value of α. Of course, the system’s coercivity gets weaker and weaker when the parameter α is getting closer to zero from positive. For small values of α, the continuity of the velocity approximation is less enforced. This can be seen from the numerical solution in Figures 5.5 and 5.6. But the discontinuity is suppressed when the value of α gets large as shown in Figure 5.7. • In the streamline portrait, the primary eddy and two corner eddies are clearly visible. • The non-symmetric finite element formulation is quite stable with respect to the parameter value of α. But extra care might be needed for solving 16
Fig. 5.3. The first velocity component profile for the lid driven cavity problem, obtained by using the symmetric formulation with various values of α. velocity u1, alpha=8
velocity u1, alpha=128
1
1
0.5
0.5
0
0
−0.5 1
−0.5 1
0 0.5
0 0.5
0.5 0
0.5 0
1
1
velocity u1, alpha=512
1 0.5 0 −0.5 1
0 0.5
0.5 0
1
the resulting non-symmetric matrix problem. GMRES was employed in our numerical experiments and the method has an acceptable convergence for the size of problem illustrated in this article. • The condition number for the resulting matrix is proportional to h−4 where h is the mesh size of the finite element partition. Therefore, neither CG nor GMRES provides a fast solver for the matrix problem. A good preconditioner is clearly needed in order to speedup the performance of CG and GMRES. Preconditioning issues will be further explored in forthcoming papers.
Fig. 5.6. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 0.1. velocity u1, alpha=0.1 (NS)
velocity u2, alpha=0.1 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5 0
0 0.5
0.5
0.5 0
1
17
1
Fig. 5.4. The second velocity component profile for the lid driven cavity problem, obtained by using the symmetric formulation with various values of α. velocity u2, alpha=8
velocity u2, alpha=128
0.6
0.6
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.4
−0.4
1
1
0 0.5
0 0.5
0.5 0
0.5 0
1
1
velocity u2, alpha=512
0.6 0.4 0.2 0 −0.2 −0.4 1
0 0.5
0.5 0
1
Fig. 5.5. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 0.01. velocity u1, alpha=0.01 (NS)
velocity u2, alpha=0.01 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5 0
0 0.5
0.5
0.5 0
1
1
Fig. 5.7. The velocity profile for the lid driven cavity problem, obtained by using the nonsymmetric formulation with α = 1. velocity u1, alpha=1 (NS)
velocity u2, alpha=1 (NS)
0.6 1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.5 1
1
0 0.5
0.5 0
1
0 0.5
0.5
18 0
1
REFERENCES [1] R.A. Adams, J.J.F. Fournier, Sobolev Spaces, Elsevier, 2nd ed., Oxford, 2003. [2] Douglas N. Arnold, Richard S. Falk and Ragnar Winther, Multigrid in H(div) and H(curl), Numer. Math., 85 (2000) 197–218. [3] F. Brezzi, J. Douglas, and L. Marini, Two families of mixed finite elements for second order elliptic problem, Numer. Math., 47 (1985) 217–235. [4] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Elements, Springer-Verlag, New York, 1991. [5] Z. Cai, R. Parashkevov, T. Russell and X. Ye, Domain decomposition for a mixed finite element method in three dimensions, SIAM J. Numer. Anal., 41 (2003) 181–194. [6] B. Cockburn, G. Kanschat, and D. Schotzau, A locally conservative LDG method for the incompressible Navier-Stokes equations, Math. Comp., 74 (2005) 1067–1095. ¨ tzau, and C. Schwab, Local discontinuous Galerkin [7] B. Cockburn, G. Kanschat, D. Scho methods for the Stokes system, SIAM J. Numer. Anal., 40 (2002) 319–343. [8] M. Costabel and M. Dauge, Crack singularities for general elliptic systems, Math. Nachr., 235 (2002) 29–49. [9] M. Dauge, Elliptic Boundary Value Problems in Corner Domains - Smoothness and Asymptotics of Solutions, Lecture Notes in Math. 1341, Springer-Verlag, Berlin, 1988. [10] M. Dauge, Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. Part I. Linearized Equations, SIAM J. Math. Anal., 20 (1989) 74–97. [11] J. Douglas, Z. Cai and X. Ye, A Stable Quadrilateral nonconforming Element for the NavierStokes Equations, Calcolo, 36 (1999) 215–232. [12] R.E. Ewing and J. Wang, Analysis of the Schwarz algorithm for mixed finite element methods, R.A.I.R.O. Modelisation Mathematique Analyse Numerique, 26 (1992) 739–756. [13] R.E. Ewing and J. Wang, Analysis of multilevel decomposition iterative methods for mixed finite element methods, R.A.I.R.O. Mathematical Modeling and Numerical Analysis, 28 (1994) 377–398. [14] V. Girault and J.-L. Lions, Two-grid finite-element schemes for the steady Navier-Stokes problem in polyhedra, Port. Math. (N.S.), 58 (2001) 25–57. [15] V. Girault and P.A. Raviart, Finite Element Methods for the Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, 1986. `re, and M.F. Wheeler, A discontinuous Gelerkin method with noncon[16] V. Girault, B. Rivie forming domain decomposition for Stokes and Navier-Stokes problems, Math. Comp., 74 (2004) 53–84. [17] D. Griffiths, Finite element for incompressible flow, Math. Meth. in Appl. Sci., 1 (1979), 16–31. [18] D. Griffiths, The construction of approximately divergence-free finite element, The Mathematics of Finite Element an Its Applications III, Ed. J.R. Whiteman, Academic Press, 1979. [19] D. Griffiths, An approximately divergence-free 9-node velocity element for incompressible flows, Inter. J. Num. Meth. in Fluid, 1 (1981) 323–346. [20] P. Grisvard, Boundary Value Problems in Non-Smooth Domains, Pitman, London, 1985. [21] K. Gustafson and R. Hartman, Divergence-free basis for finite element schemes in hydrodynamics, SIAM J. Numer. Anal., 20 (1983) 697–721. [22] K. Gustafson and R. Hartman, Graph theory and fluid dynamics, SIAM J. Alg. Disc. Meth., 6 (1985) 643–656. [23] M. D. Gunzburger, Finite Element Methods for Viscous Incompressible Flows, A Guide to Theory, Practice and Algorithms, Academic, San Diego, 1989. [24] O. A. Karakashian and W. N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM J. Numerical Analysis, 35 (1998) 93–120. [25] J. R. Kweon, Corner singularity for compressible Navier-Stokes problem, in Proceedings of the International Congress of Mathematics, Higher Education Press, China, 2002. [26] J. Liu and C. Shu, A high order discontinuous Galerkin method for 2D incompressible flows, J. Comput. Phys., 160 (2000) 577–596. [27] P.A. Raviart and J.M. Thomas, A mixed finite element method for second order elliptic problems, Mathematical aspects of the finite element method, Eds. I. Galligani and E. Magenes, Lecture notes in Mathematics, Vol. 606., Springer-Verlag, 1977. [28] R. Scheichl, Decoupling three-dimensional mixed problems using divergence-free finite elements, SIAM J. Sci. Comp., 23 (2002) 1752–1776. [29] P.N. Shankar and M.D. Deshpande, Fluid Mechanics in the driven cavity, Annu. Rev. Fluid Mech., 32 (2000) 93–136. 19
[30] J. Wang and X. Ye, New finite element methods in computational fluid dynamics by H(div) elements, SIAM J. Numer. Anal., 45 (2007) 1269–1286. [31] J. Wang, X. Wang and X. Ye, Finite Element Methods for the Navier-Stokes equations by H(div) Elements, Journal of Computational Mathematics, 26 (2008) 410–436. [32] X. Ye and C. Hall, A Discrete divergence free basis for finite element methods, Numerical Algorithms, 16 (1997) 365–380. [33] X. Ye and C. Hall, The Construction of null basis for a discrete divergence operator, J. Computational and Applied Mathematics, 58 (1995) 117–133. [34] X. Ye, Discontinuous stable elements for the incompressible flow, Advances in Computational Mathematics, 20 (2004) 333–345. [35] X. Ye, Construction divergence free space for incompressible Navier-Stokes equations, Ph.D thesis, 1990.
20
