- Email: [email protected]
Dual-mixed finite elements for the three-field Stokes model as a finite volume method on staggered grids
Available online at www.sciencedirect.com
ScienceDirect
This space is reserved for the Procedia header, do not use it This space reserved for the header, do not use it ProcediaisComputer Science 108CProcedia (2017) 2265–2274 This space is reserved for the Procedia header, do not use it
International Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland
Dual-mixed finite elements for the three-field Stokes model Dual-mixed finite elements for the three-field Stokes model as a finite method staggeredStokes grids Dual-mixed finitevolume elements for theon as a finite volume method onthree-field staggered grids model 1 Jisheng Kou and Shuyu Sun2 as a finite volume method on staggered grids 1 Jisheng Kou and Shuyu Sun2
1
School of Mathematics and Statistics, Hubei Engineering University, 1 2 Xiaogan 432000, Hubei, China. JishengHubei Kou and Division Shuyu Sun Computational Transport Phenomena Laboratory, of Physical Science and Engineering, School of Mathematics and Statistics, Engineering University, Xiaogan 432000, Hubei, China. 2 Computational Transport Phenomena Laboratory, Thuwal Division23955-6900, of Physical Kingdom Science and Engineering, King Abdullah University of Science and Technology, of Saudi Arabia. 1 School of Mathematics Statistics, Hubei Engineering University, Xiaogan 432000, Hubei,Arabia. China. [email protected] King Abdullah Universityand of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi 2 Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, [email protected] King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia. [email protected] 1 2
Abstract Abstract In this paper, a new three-field weak formulation for Stokes problems is developed, and from this, In this paper, afinite new element three-field weak formulation problems is developed, and from this, a dual-mixed method is proposed for on Stokes a rectangular mesh. In the proposed mixed Abstract a dual-mixed elementofmethod proposed on a rectangular In the proposed mixed methods, the finite components stress is tensor are approximated by mesh. piecewise constant functions In this paper, a new three-field weak formulation for Stokes problems is developed, and from this, methods, the components stress and tensor are approximated by by piecewise constant functions or Q1 functions, while the of velocity pressure are discretized the lowest-order Raviarta dual-mixed finite element method is proposed on a rectangular mesh. In the proposed mixed or Q1 functions, while velocityconstant and pressure are discretized byUsing the lowest-order RaviartThomas element and thethe piecewise functions, respectively. quadrature rules, we methods, the components of stress tensor are approximated by piecewise constant functions Thomas element piecewise functions, respectively. Using quadrature rules, we demonstrate thatand thisthe scheme can constant be reduced into a finite volume method on staggered grid, or Q1 functions, while the velocity and pressure are discretized by the lowest-order Raviartdemonstrate that this scheme can be reduced a finite and volume method on staggered grid, which is extensively used in computational fluidinto mechanics engineering. Thomas element and the piecewise constant functions, respectively. Using quadrature rules, we which is extensively used in computational fluid mechanics and engineering. Keywords: Stokes three-field model, mixed elements, staggered demonstrate thatequations, this scheme can beB.V. reduced intofinite a finite volume methodgrid. on staggered grid, © 2017 The Authors. Published by Elsevier Keywords: Stokes equations, model, mixed finite elements, staggeredon grid. Peer-review under responsibility ofthree-field the scientific committee of the International Computational Science which is extensively used in computational fluid mechanics and Conference engineering. Keywords: Stokes equations, three-field model, mixed finite elements, staggered grid.
1 Introduction 1 Introduction It is well-known that Stokes flow is a typical fluid flow that is dominated by viscous forces. 1 Introduction It is well-known Stokes flow is a typical fluid flowof that is dominated Here, we use the that notation η to represent the viscosity a fluid. Denote a by 2D viscous domain forces. by Ω. T notation η to represent the viscosity of a fluid. Denote a 2D domain by Ω. Here,u we usev]the Let = [u, and p denote the flow velocity and pressure, respectively. We consider the It is well-known that Stokes flow is a typical fluid flow that is dominated by viscous forces. Let u = [u, v]T steady and p Stokes denoteequation the flowasvelocity and pressure, respectively. We consider the incompressible, Here, we use the notation η to represent the viscosity of a fluid. Denote a 2D domain by Ω. incompressible,T steady Stokes equation as Let u = [u, v] and p denote−∇ the· (η∇u) flow velocity pressure, + ∇p =and f, x ∈ Ω, respectively. We consider the (1) incompressible, steady Stokes−∇ equation as+ ∇p = f , x ∈ Ω, · (η∇u) (1) along with the incompressibility condition (i.e. mass balance equation) −∇condition · (η∇u) +(i.e. ∇p = f , balance x ∈ Ω, equation) (1) along with the incompressibility mass ∇ · u = 0, x ∈ Ω. (2) along with the incompressibility condition mass ∇ · u =(i.e. 0, x ∈ Ω.balance equation) (2) Here, x = [x, y]T , and f = [fx , fy ]T is the source term. The boundary condition is given as T∇ · u = 0, x ∈ Ω. Here, y]T For , andsake f =of[fsimplicity, source term. condition is given as u = uxB = on[x, ∂Ω. we assume uB = The 0; inboundary fact, for u have (2) its x , fy ] is the B �= 0, we u = uB onextension ∂Ω.T For on sake simplicity, wedenoted assume by uBuB= yet, 0; in fact, for case, uB �=we 0,just we consider have its continuous theofentire domain, and in this T Here, x = [x, y] , and f = [fx , fy ] is the source term. The boundary condition is given as continuous domain, denoted by u ayet, and in this case, wefor just consider ˆ = u − uBextension ˆ on u , where u = the 0 onentire the boundary. In addition, compatible condition pressure is u = uB on ∂Ω. For sake of simplicity, we assume uB B= 0; in fact, for uB �= 0, we have its ˆ = u − ufor ˆ = Ω0p(x)dx u where u on the = boundary. In addition, a compatible condition for pressure is required, 0. B , example, continuous extension on the entire domain, denoted by uB yet, and in this case, we just consider required, for example, Ω p(x)dx = 0. ˆ = u − uB , where u ˆ = 0 on the boundary. In addition, a compatible condition for pressure is u 1 required, for example, Ω p(x)dx = 0. 1
1877-0509 © 2017 The Authors. Published by Elsevier B.V. Peer-review under responsibility of the scientific committee of the International Conference on Computational Science 10.1016/j.procs.2017.05.093
1
2266
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
In the formulation (1), the velocity and pressure are treated as the two primal variables, so it is called two-field model. The three-field model of Stokes flow introduces the stress tensor as a new variable except for velocity and pressure. We define the stress tensor τxx τxy τ = . τyx τyy Then (1) can be reformulated as −∇ · τ + ∇p = f ,
η −1 τ = ∇u.
(3)
The three-field model has been extensively used in nonlinear Stokes problems [2, 6, 12]. The mixed finite element method is a very popular method for simulating flow problems in porous media [1,9], Stokes and Navier-Stokes equations [3,4,7]. In particular, a series of research activities is carried out on simple and efficient nonconforming finite element methods for the two-field model of Stokes problems, for example, [5, 10, 11]. These mixed methods discretize the velocity and pressure in different finite element spaces. On the other hand, the dual-mixed finite element methods for a three-field model of Stokes problems have also been developed, for example, [2,6,12]. Compared to mixed methods based on two-field model, dual-mixed methods for a three-field model provide more flexibility in approximating three variables using different finite element spaces, so they may achieve higher computational accuracy for not only the pair of velocity and pressure but also the stress. The relationship between mixed finite element methods and finite difference or finite volume schemes is an interesting research issue. It has been shown in [1] that the cell-centered finite difference method for elliptic problems can be derived from the lowest-order Raviart-Thomas (RT0 ) spaces [13] on rectangles using appropriate quadrature rules. The staggered grid introduced by Harlow and Welch [8] has been used extensively and popularly in computational fluid mechanics. Finite volume schemes on staggered grid have also been the basis of many industrial applications. It is shown in [7] that the MAC method can be derived by combining a velocity-vorticity mixed finite element method with an adequate quadrature formula. In this paper, we will first construct a three-field weak formulation for Stokes problems, and from this, a new dual-mixed finite element method is proposed. We discretize the velocity in the lowest-order Raviart-Thomas space and the pressure with piecewise constant functions, respectively. The components of stress tensor are approximated by piecewise constant functions or Q1 functions. Applying trapezoid quadrature rules to our proposed scheme, we derive a finite volume method on staggered grid, which is extensively used in computational fluid mechanics and engineering.
2
Finite volume method on staggered grid
We now formulate a finite volume method on staggered grid, which is one of the most common methods for simulating Stokes equations in practical applications. For sake of simplification, we focus on the case of rectangular domain, but the proposed scheme can be extended to the other regular domains. Moreover, the viscosity is assumed to be a constant. We consider a rectangular domain Ω = [0, Lx ] × [0, Ly ]. For sake of simplification, we use a uniform rectangular mesh for the domain, and denote the mesh size by h. The interval [0, Lx ] is uniformly divided into M subintervals, while the interval [0, Ly ] is uniformly divided into N subintervals. The variables located on the vertexes of each element are denoted by integer subscripts or superscripts, while variables on the middle points of an edge or element are denoted 2
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
by fraction subscripts or superscripts. For example, for a control volume [xi , xi+1 ] × [yj , yj+1 ], the pressure is located at its center (denoted by pi+ 12 ,j+ 12 ), the velocity component u is located at the middle of the left and right edges (denoted by ui,j+ 12 and ui+1,j+ 12 ), and the velocity components v is located at the middle of the top and bottom edges (denoted by vi+ 21 ,j and 1 1 1 1 i+ 1 ,j+ 1 T i+ 1 ,j+ 1 vi+ 21 ,j+1 ). We define the discrete source term f i+ 2 ,j+ 2 = fx 2 2 , fy 2 2 as f i+ 2 ,j+ 2 = xi+1 yj+1 1 f (x, y)dydx. The discrete form of Stokes equation on the interior elements is h2 x i yj expressed as η 4ui,j+ 12 − ui−1,j+ 12 − ui+1,j+ 21 − ui,j− 21 − ui,j+ 32 1 i− 1 ,j+ 1 i+ 1 ,j+ 1 (4) + h pi+ 21 ,j+ 12 − pi− 12 ,j+ 12 = h2 fx 2 2 + fx 2 2 , 2 η 4vi+ 21 ,j − vi− 12 ,j − vi+ 32 ,j − vi+ 21 ,j−1 − vi+ 21 ,j+1 1 i+ 1 ,j+ 1 i+ 1 ,j− 1 + h pi+ 12 ,j+ 12 − pi+ 12 ,j− 12 = h2 fy 2 2 + fy 2 2 . 2
(5)
The mass balance equation on the interior elements is discretized as ui+1,j+ 12 − ui,j+ 21 + vi+ 12 ,j+1 − vi+ 21 ,j = 0.
(6)
The discrete forms on the boundary elements can be similarly obtained by substituting the boundary conditions.
3
Dual-mixed finite elements for Stokes equation
We define a function space as 2 2×2 ∂ τˆxy ∂ τˆyx τˆxx τˆxy 2 2 ∈ L (Ω), ∈ L (Ω) . ∈ L (Ω) : S = τˆ = τˆyx τˆyy ∂x ∂y A bilinear operator is defined as ∂τyx ∂v ∂τxy ∂u τxx − v −u + τyy dx. A∇ (u, τ ) = ∂x ∂x ∂y ∂y Ω
(7)
A mixed weak form of Stokes equation is to seek (τ , p, u) ∈ S × L2 (Ω) × H(div; Ω) such that A∇ (v, τ ) − (p, ∇ · v) = (f , v), v ∈ H(div; Ω),
(8)
η −1 τ , τˆ = A∇ (u, τˆ ), τˆ ∈ S,
(9)
(∇ · u, w) = 0, w ∈ L2 (Ω).
(10)
We use a uniform rectangular mesh, which is the same to that stated in the above section. Let K be any element. The finite element space S h is defined as τ : τˆxx ∈ P0 (K) , τˆyy ∈ P0 (K) , τˆxy ∈ Q1 (K) , τˆyx ∈ Q1 (K)}. S h = {ˆ
(11) 3
2267
2268
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
The finite element space RT0 on an element K of a rectangular mesh is defined as � � RT0 = {v : v = a1K + a2K x, a3K + a4K y , aiK ∈ R, 1 ≤ i ≤ 4}.
(12)
More precisely, the basis function of RT0 can be expressed as v = [φxi,j+ 1 (x, y), 0]T or v = 2
[0, φyi+ 1 ,j (x, y)]T , where 2
φxi,j+ 1 (x, y) = 2
xi+1 −x xi+1 −xi , x−xi−1 xi −xi−1 ,
yj+1 −y yj+1 −yj , y−yj−1 yj −yj−1 ,
φyi+ 1 ,j (x, y) 2
=
0,
0,
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi−1 , xi ], y ∈ [yj , yj+1 ] , elsewhere
(13)
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi , xi+1 ], y ∈ [yj−1 , yj ] . elsewhere
(14)
Thus, the velocity uh ∈ RT0 has the unique representation uh (x) =
N M+1 ��
ui,j+ 12 φxi,j+ 1 (x), vh (x) = 2
i=1 j=1
+1 M N � �
vi+ 12 ,j φyi+ 1 ,j (x), x ∈ Ω.
The basis functions of the finite element space S h have the following forms � 1, x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] i+ 1 ,j+ 1 i+ 1 ,j+ 1 , τˆxx 2 2 (x, y) = τˆyy 2 2 (x, y) = 0, elsewhere
i,j i,j (x, y) = τˆxy (x, y) = τˆyx
xi+1 −x xi+1 −xi x−xi−1 xi −xi−1 xi+1 −x xi+1 −xi x−xi−1 xi −xi−1
yj+1 −y yj+1 −yj , yj+1 −y yj+1 −yj , y−yj−1 yj −yj−1 , y−yj−1 yj −yj−1 ,
0,
(15)
2
i=1 j=1
(16)
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi−1 , xi ], y ∈ [yj , yj+1 ] x ∈ [xi , xi+1 ], y ∈ [yj−1 , yj ] , x ∈ [xi−1 , xi ], y ∈ [yj−1 , yj ] elsewhere
(17)
We define the finite element approximation of the stress as � � τxx,h τxy,h , x ∈ Ω, τ h (x) = τyx,h τyy,h
(18)
where τxx,h (x) =
N M � �
i+ 1 ,j+ 12
τxx 2
i+ 1 ,j+ 12
τˆxx 2
(x),
τyy,h (x) =
i=1 j=1
τxy,h (x) =
N M � �
i+ 1 ,j+ 12
τˆyy 2
(x),
(19)
i=1 j=1
+1 M+1 � N� i=1 j=1
i,j i,j τxy τˆxy (x),
τyx,h (x) =
+1 M+1 �N � i=1 j=1
We define the finite element space Wh on an element K as Wh = {w : w is constant in K}. 4
i+ 1 ,j+ 12
τyy 2
i,j i,j τyx τˆyx (x).
(20)
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
Furthermore, we define the basis functions of Wh as 1, x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] . wi+ 21 ,j+ 12 (x) = 0, elsewhere
(21)
The discrete forms of the pressure and source term are expressed as ph (x) =
N M
pi+ 21 ,j+ 21 wi+ 12 ,j+ 21 (x), f h (x) =
N M
1
1
f i+ 2 ,j+ 2 wi+ 12 ,j+ 12 (x), x ∈ Ω. (22)
i=1 j=1
i=1 j=1
The mixed finite element formulation is to seek (τ h , ph , uh ) ∈ S h × Wh × RT0 such that A∇ (v, τ h ) − (ph , ∇ · v) = (f h , v), v ∈ RT0 ,
(23)
−1 η τ h , τˆ = A∇ (uh , τˆ ), τˆ ∈ S h ,
(24)
(∇ · uh , w) = 0, w ∈ Wh .
(25)
From the properties of mixed finite elements, the proposed scheme has approximate errors of at least one order, but because of limitation in the length of this article, the detailed analysis will be presented in our future work.
4
Equivalence forms
We now demonstrate the equivalence between the proposed finite element method with appropriate quadrature rules and the finite volume method on staggered grid stated in Section 2. For simplification, we derive this equivalence for the interior elements only. i+ 21 ,j+ 21 −1 0 τ ˆ xx , we Firstly, we consider the representation of η τ h , τˆ . Taking τˆ = 0 0 have −1 i+ 1 ,j+ 1 η τ h , τˆ = h2 η −1 τxx 2 2 . (26) 0 0 leads to The choice of τˆ = i+ 1 ,j+ 1 0 τˆyy 2 2
Let τˆ =
i,j 0 τˆxy 0 0
η −1 τ h , τˆ
−1 i+ 1 ,j+ 1 η τ h , τˆ = h2 η −1 τyy 2 2 .
(27)
, we obtain
−1 i,j η τxy,h , τˆxy xi+1 yj+1 i,j = η −1 τxy,h τˆxy dydx + =
xi
+
yj
xi+1
xi
yj
i,j η −1 τxy,h τˆxy dydx + yj−1
xi
xi−1 xi
xi−1
yj+1 i,j η −1 τxy,h τˆxy dydx
yj
yj i,j η −1 τxy,h τˆxy dydx.
(28)
yj−1
5
2269
2270
Dual-mixed finite elements for Jisheng the three-field model Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274
Kou and Sun
We apply trapezoid quadrature rule to the first term on the right-hand side of (28) � xi+1 � yj+1 i,j η −1 τxy,h τˆxy dydx xi
=
�
yj
xi+1 �
xi
=
≈
yj+1
yj
j+1 i+1 � � k,l k,l i,j η −1 τxy τˆxy τˆxy dydx k=i l=j
�2 yj+1 − y dx dy yj+1 − yj xi yj �2 � xi+1 � yj+1 � xi+1 − x x − xi yj+1 − y −1 i+1,j dx dy +η τxy xi+1 − xi xi+1 − xi yj+1 − yj xi yj �2 � yj+1 � xi+1 � yj+1 − y y − yj xi+1 − x i,j+1 dx dy +η −1 τxy x − x y i+1 i j+1 − yj yj+1 − yj yj xi � xi+1 � yj+1 yj+1 − y y − yj xi+1 − x x − xi −1 i+1,j+1 dx dy +η τxy x − x x − x y i+1 i i+1 i j+1 − yj yj+1 − yj yj xi i,j η −1 τxy
�
xi+1
�
xi+1 − x xi+1 − xi
�2
�
yj+1
�
h2 i,j τ . 4η xy
(29)
The rest three terms on the right-hand side of (28) can be approximated by the similar way, and consequently, in this case, we get � −1 � � � i,j i,j η τ h , τˆ ≈ h2 η −1 τxy = η −1 τxy,h , τˆxy . (30) By the similar appraoches, we take τˆ =
�
0 i,j τˆyx
0 0
�
and obtian
� � � i,j i,j η −1 τ h , τˆ = η −1 τyx,h , τˆyx ≈ h2 η −1 τyx .
�
Secondly, we consider the form of A∇ (uh , τˆ ). We take � � i+ 1 ,j+ 1 � τˆxx τˆxy τˆxx 2 2 = τˆ = τˆyx τˆyy 0
0 0
�
(31)
,
and then A∇ (uh , τˆ ) becomes A∇ (uh , τˆ )
Similarly, let τˆ =
�
� � � xi+1 � yj+1 � ∂uh ∂uh i+ 12 ,j+ 21 τˆxx dx = τˆxx dx ∂x ∂x Ω xi yj � xi+1 � yj+1 � � ∂uh dx = h ui+1,j+ 21 − ui,j+ 12 . = ∂x xi yj
=
0
� �
�
0 i+ 21 ,j+ 12
0 τˆyy
A∇ (uh , τˆ ) =
�
6
, we get
xi+1
xi
(32)
�
yj+1
yj
� � ∂vh dx = h vi+ 12 ,j+1 − vi+ 12 ,j . ∂x
(33)
Dual-mixed finite elements for Jisheng the three-field model Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274
Taking τˆ =
i,j 0 τˆxy 0 0
Kou and Sun
gives
i,j ∂ τˆxy dx −vh ∂x Ω xi+1 yj+1 xi yj+1 i,j i,j ∂ τˆxy ∂ τˆxy − dydx − dydx vh vh ∂x ∂x xi yj xi−1 yj xi yj xi+1 yj i,j i,j ∂ τˆxy ∂ τˆxy dydx − dydx − vh vh ∂x ∂x xi yj−1 xi−1 yj−1 1 yj+1 yj+1 − y y − yj (yj+1 − y) dy vi+ 21 ,j + vi+ 21 ,j+1 h yj yj+1 − yj yj+1 − yj 1 yj+1 yj+1 − y y − yj vi− 21 ,j (yj+1 − y) dy − + vi− 21 ,j+1 h yj yj+1 − yj yj+1 − yj 1 yj yj − y y − yj−1 vi+ 21 ,j−1 (y − yj−1 ) dy + + vi+ 21 ,j h yj−1 yj − yj−1 yj − yj−1 1 yj yj − y y − yj−1 (y − yj−1 ) dy. v 1 − + vi− 21 ,j h yj−1 i− 2 ,j−1 yj − yj−1 yj − yj−1
A∇ (uh , τˆ ) = =
=
(34)
Applying trapezoid quadrature rule to the terms on the right-hand side of (34), we deduce A∇ (uh , τˆ ) ≈ h vi+ 21 ,j − vi− 12 ,j . (35) Taking τˆ = deduce
0 i,j τˆyx
0 0
, and furthermore, by the similar approaches used in (34)-(35), we A∇ (uh , τˆ ) ≈ h ui,j+ 21 − ui,j− 12 .
(36)
Combining the above derivations, we derive the relations between τ h and uh as η η i+ 1 ,j+ 1 i+ 1 ,j+ 1 ui+1,j+ 12 − ui,j+ 12 , τyy 2 2 = vi+ 12 ,j+1 − vi+ 21 ,j , τxx 2 2 = h h η ui,j+ 12 − ui,j− 21 . h T Thirdly, we consider the expressions of A∇ (v, τ h ). Take v = φxi,j+ 1 , 0 i,j ≈ τxy
η vi+ 12 ,j − vi− 21 ,j , h
i,j τyx ≈
(37)
(38)
2
A∇ (v, τ h )
=
∂φx 1 i,j+ 2
Ω
∂x
τxx,h −
φxi,j+ 1 2
∂τyx,h ∂y
dx.
(39)
Using their definitions, the first term on the right-hand side of (39) is calculated as ∂φx 1 xi+1 yj+1 ∂φx 1 i,j+ 2 i+ 12 ,j+ 12 i,j+ 2 τxx,h dx = τxx dydx ∂x ∂x xi yj Ω 7
2271
2272
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model
+
xi
xi−1
=
yj+1
yj
∂φxi,j+ 1
i− 1 ,j+ 1 τxx 2 2
2
∂x
Kou and Sun
dydx
i+ 1 ,j+ 1 i− 1 ,j+ 1 −h τxx 2 2 − τxx 2 2 .
(40)
The second term on the right-hand side of (39) is rewritten as xi+1 yj+1 xi yj+1 ∂τyx,h ∂τyx,h ∂τyx,h dx = dydx + dydx φxi,j+ 1 φxi,j+ 1 φxi,j+ 1 2 2 2 ∂y ∂y ∂y Ω xi yj xi−1 yj xi+1 i,j+1 xi+1 − x xi+1 − x i,j = dx τyx − τyx xi+1 − xi xi+1 − xi xi xi+1 i+1,j+1 xi+1 − x x − xi i+1,j dx − τyx τyx + xi+1 − xi xi+1 − xi x xi i x − xi−1 i,j+1 x − xi−1 i,j dx τyx − τyx + x − x x i i−1 i − xi−1 xi−1 xi i−1,j+1 x − xi−1 xi − x i−1,j + dx. (41) τyx − τyx xi − xi−1 xi − xi−1 xi−1 Applying trapezoid quadrature rule to the four terms on the right-hand side of (41) leads to i,j+1 ∂τyx,h i,j dx ≈ −h τyx −φxi,j+ 1 . (42) − τyx 2 ∂y Ω T Let v = 0, φyi+ 1 ,j , we have 2
A∇ (v, τ h ) =
∂φy 1 i+ ,j 2
∂y
Ω
τyy,h −
φyi+ 1 ,j 2
∂τxy,h ∂x
(43)
dx.
The first term on the right-hand side of (43) is directly calculated as ∂φy 1 i+ ,j 2
Ω
∂y
τyy,h dx =
xi+1
yj
xi+1
xi
+
2
yj−1
xi
=
∂φyi+ 1 ,j yj+1
yj
∂y
i+ 1 ,j− 12
τyy 2
∂φyi+ 1 ,j 2
∂y
dydx
i+ 1 ,j+ 12
τyy 2
i+ 1 ,j+ 1 i+ 1 ,j− 1 −h τyy 2 2 − τyy 2 2 .
dydx (44)
For the second term on the right-hand side of (43), we use the trapezoid quadrature rule xi+1 yj xi+1 yj+1 ∂τxy,h ∂τxy,h ∂τxy,h y y − dx = − dydx − dydx φi+ 1 ,j φyi+ 1 ,j φi+ 1 ,j 2 2 2 ∂x ∂x ∂x xi yj−1 xi yj Ω yj+1 i+1,j yj+1 − y yj+1 − y i,j = − dy τxy − τxy yj+1 − yj yj+1 − yj yj yj+1 i+1,j+1 yj+1 − y y − yj i,j+1 − dy τxy − τxy yj+1 − yj yj+1 − yj yj 8
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model yj
Kou and Sun
y − yj−1 y − yj−1 dy yj − yj−1 yj − yj−1 yj−1 yj i+1,j−1 y − yj−1 yj − y i,j−1 − dy τxy − τxy yj − yj−1 yj − yj−1 yj−1 i+1,j i,j ≈ −h τxy − τxy . (45) T Finally, we consider the formulations of (ph , ∇ · v) and (f h , v). Taking v = φxi,j+ 1 , 0 2 T y and v = 0, φi+ 1 ,j respectively, we obtain 2 ∂φxi,j+ 1 2 = −h pi+ 12 ,j+ 12 − pi− 21 ,j+ 12 , (46) (ph , ∇ · v) = ph , ∂x −
(ph , ∇ · v) =
ph ,
i,j i+1,j − τxy τxy
∂φyi+ 1 ,j 2
∂y
= −h pi+ 12 ,j+ 12 − pi+ 12 ,j− 12 .
(47)
T T respectively, we get For (f h , v), taking v = φxi,j+ 1 , 0 and v = 0, φyi+ 1 ,j 2
2
1 i+ 1 ,j+ 1 i− 1 ,j+ 1 (f h , v) = fx,h , φxi,j+ 1 = h2 fx 2 2 + fx 2 2 , 2 2
(48)
1 i+ 1 ,j+ 1 i+ 1 ,j− 1 (f h , v) = fy,h φyi+ 1 ,j = h2 fy 2 2 + fy 2 2 . 2 2
(49)
Thus the equation (24) can be rewritten as i− 1 ,j+ 1 i+ 1 ,j+ 1 i,j+1 i,j −h τxx 2 2 − τxx 2 2 + τyx + h pi+ 21 ,j+ 12 − pi− 12 ,j+ 21 − τyx 1 i+ 1 ,j+ 1 i− 1 ,j+ 1 = h2 f x 2 2 + f x 2 2 , 2 i+ 1 ,j+ 1 i+ 1 ,j− 1 i+1,j i,j −h τyy 2 2 − τyy 2 2 + τxy + h pi+ 21 ,j+ 12 − pi+ 12 ,j− 12 − τxy 1 i+ 1 ,j+ 1 i+ 1 ,j− 1 = h2 f y 2 2 + f y 2 2 . 2
(50)
(51)
The discrete forms (4) and (5) are obtained substituting (37) and (38) into the equations (50) and (51). For the equation (∇ · uh , w) = 0, let w = 1 on an element K, we get (6). We note that the above derivations are carried out on the interior elements, but the boundary elements can be handled by the similar approaches and taking into account the boundary conditions.
5
Numerical results
We test the benchmark driven cavity problem. For Ω = [0, 1]2 , we have that uB = [1, 0]T on the top boundary (that is, the portion of the boundary satisfying 0 ≤ x ≤ 1 and y = 1) and uB = 0 on the rest portions of the boundary. The source term is zero. The viscosity is set as η = 1. Figs. 1 shows the computed flow quiver on a rectangular mesh composed of 30 × 30 square elements. 9
2273
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model
Kou and Sun
Flow quiver
1 0.8 0.6
Y
2274
0.4 0.2 0 0
0.2
0.4
X
0.6
0.8
1
Figure 1: The flow quiver.
6
Conclusions
A new weak formulation is developed for the three-field model of Stokes problems. Based on this formulation, we proposed a dual-mixed finite element method on a rectangular mesh. A finite volume method on staggered grid is derived by applying trapezoid quadrature rules for the proposed scheme.
References [1] T. Arbogast, M.F. Wheeler and I. Yotov. Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences. SIAM J. Numer. Anal., 34: 828–852, 1997. [2] J. Baranger, K. Najib and D. Sandri. Numerical analysis of a three-beads model for a quasiNewtonian flow. Comput. Methods Appl. Mech. Engrg., 109: 281–292, 1993. [3] S. C. Brenner and R. Scott. The mathematical theory of finite element methods, volume 15. Springer, 2007. [4] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer Verlag, 1991. [5] M. Crouzeix and P.A. Raviart. Conforming and nonconforming finite element methods for solving the stationary stokes equations. RAIRO Analysis Numerical, 7:33–36, 1973. [6] V. J. Ervin, J. S. Howell and I. Stanculescu. A dual-mixed approximation method for a threefield model of a nonlinear generalized Stokes problem. Comput. Methods Appl. Mech. Engrg., 197: 2886–2900, 2008. [7] V. Girault, H. Lopez. Finite-element error estimates for the MAC scheme. IMA Journal of Numerical Analysis, 16(3): 347-379, 1996. [8] F. H. Harlow and J. E. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys. Fluids, 8: 2182–2189, 1965. [9] J. Kou, S. Sun and Y. Wu. Mixed finite element-based fully conservative methods for simulating wormhole propagation. Comput. Methods Appl. Mech. Engrg., 298: 279–302, 2016. [10] B. P. Lamichhane. A nonconforming finite element method for the Stokes equations using the Crouzeix-Raviart element for the velocity and the standard linear element for the pressure. Int. J. Numer. Meth. Fluids, 74: 222–228, 2014. [11] J. Li and Z. Chen. A new local stabilized nonconforming finite element method for the Stokes equations. Computing, 82:157–170, 2008. [12] H. Manouzi and M. Farhloul. Mixed finite element analysis of a nonlinear three-fields Stokes model. IMA J. Numer. Anal., 21: 143–164, 2001. [13] R.A. Raviart and J.M. Thomas. Primal hybrid finite element methods for 2nd order elliptic equations. Math. Comput., 31: 391–413, 1977.
10
ScienceDirect
This space is reserved for the Procedia header, do not use it This space reserved for the header, do not use it ProcediaisComputer Science 108CProcedia (2017) 2265–2274 This space is reserved for the Procedia header, do not use it
International Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland
Dual-mixed finite elements for the three-field Stokes model Dual-mixed finite elements for the three-field Stokes model as a finite method staggeredStokes grids Dual-mixed finitevolume elements for theon as a finite volume method onthree-field staggered grids model 1 Jisheng Kou and Shuyu Sun2 as a finite volume method on staggered grids 1 Jisheng Kou and Shuyu Sun2
1
School of Mathematics and Statistics, Hubei Engineering University, 1 2 Xiaogan 432000, Hubei, China. JishengHubei Kou and Division Shuyu Sun Computational Transport Phenomena Laboratory, of Physical Science and Engineering, School of Mathematics and Statistics, Engineering University, Xiaogan 432000, Hubei, China. 2 Computational Transport Phenomena Laboratory, Thuwal Division23955-6900, of Physical Kingdom Science and Engineering, King Abdullah University of Science and Technology, of Saudi Arabia. 1 School of Mathematics Statistics, Hubei Engineering University, Xiaogan 432000, Hubei,Arabia. China. [email protected] King Abdullah Universityand of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi 2 Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, [email protected] King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia. [email protected] 1 2
Abstract Abstract In this paper, a new three-field weak formulation for Stokes problems is developed, and from this, In this paper, afinite new element three-field weak formulation problems is developed, and from this, a dual-mixed method is proposed for on Stokes a rectangular mesh. In the proposed mixed Abstract a dual-mixed elementofmethod proposed on a rectangular In the proposed mixed methods, the finite components stress is tensor are approximated by mesh. piecewise constant functions In this paper, a new three-field weak formulation for Stokes problems is developed, and from this, methods, the components stress and tensor are approximated by by piecewise constant functions or Q1 functions, while the of velocity pressure are discretized the lowest-order Raviarta dual-mixed finite element method is proposed on a rectangular mesh. In the proposed mixed or Q1 functions, while velocityconstant and pressure are discretized byUsing the lowest-order RaviartThomas element and thethe piecewise functions, respectively. quadrature rules, we methods, the components of stress tensor are approximated by piecewise constant functions Thomas element piecewise functions, respectively. Using quadrature rules, we demonstrate thatand thisthe scheme can constant be reduced into a finite volume method on staggered grid, or Q1 functions, while the velocity and pressure are discretized by the lowest-order Raviartdemonstrate that this scheme can be reduced a finite and volume method on staggered grid, which is extensively used in computational fluidinto mechanics engineering. Thomas element and the piecewise constant functions, respectively. Using quadrature rules, we which is extensively used in computational fluid mechanics and engineering. Keywords: Stokes three-field model, mixed elements, staggered demonstrate thatequations, this scheme can beB.V. reduced intofinite a finite volume methodgrid. on staggered grid, © 2017 The Authors. Published by Elsevier Keywords: Stokes equations, model, mixed finite elements, staggeredon grid. Peer-review under responsibility ofthree-field the scientific committee of the International Computational Science which is extensively used in computational fluid mechanics and Conference engineering. Keywords: Stokes equations, three-field model, mixed finite elements, staggered grid.
1 Introduction 1 Introduction It is well-known that Stokes flow is a typical fluid flow that is dominated by viscous forces. 1 Introduction It is well-known Stokes flow is a typical fluid flowof that is dominated Here, we use the that notation η to represent the viscosity a fluid. Denote a by 2D viscous domain forces. by Ω. T notation η to represent the viscosity of a fluid. Denote a 2D domain by Ω. Here,u we usev]the Let = [u, and p denote the flow velocity and pressure, respectively. We consider the It is well-known that Stokes flow is a typical fluid flow that is dominated by viscous forces. Let u = [u, v]T steady and p Stokes denoteequation the flowasvelocity and pressure, respectively. We consider the incompressible, Here, we use the notation η to represent the viscosity of a fluid. Denote a 2D domain by Ω. incompressible,T steady Stokes equation as Let u = [u, v] and p denote−∇ the· (η∇u) flow velocity pressure, + ∇p =and f, x ∈ Ω, respectively. We consider the (1) incompressible, steady Stokes−∇ equation as+ ∇p = f , x ∈ Ω, · (η∇u) (1) along with the incompressibility condition (i.e. mass balance equation) −∇condition · (η∇u) +(i.e. ∇p = f , balance x ∈ Ω, equation) (1) along with the incompressibility mass ∇ · u = 0, x ∈ Ω. (2) along with the incompressibility condition mass ∇ · u =(i.e. 0, x ∈ Ω.balance equation) (2) Here, x = [x, y]T , and f = [fx , fy ]T is the source term. The boundary condition is given as T∇ · u = 0, x ∈ Ω. Here, y]T For , andsake f =of[fsimplicity, source term. condition is given as u = uxB = on[x, ∂Ω. we assume uB = The 0; inboundary fact, for u have (2) its x , fy ] is the B �= 0, we u = uB onextension ∂Ω.T For on sake simplicity, wedenoted assume by uBuB= yet, 0; in fact, for case, uB �=we 0,just we consider have its continuous theofentire domain, and in this T Here, x = [x, y] , and f = [fx , fy ] is the source term. The boundary condition is given as continuous domain, denoted by u ayet, and in this case, wefor just consider ˆ = u − uBextension ˆ on u , where u = the 0 onentire the boundary. In addition, compatible condition pressure is u = uB on ∂Ω. For sake of simplicity, we assume uB B= 0; in fact, for uB �= 0, we have its ˆ = u − ufor ˆ = Ω0p(x)dx u where u on the = boundary. In addition, a compatible condition for pressure is required, 0. B , example, continuous extension on the entire domain, denoted by uB yet, and in this case, we just consider required, for example, Ω p(x)dx = 0. ˆ = u − uB , where u ˆ = 0 on the boundary. In addition, a compatible condition for pressure is u 1 required, for example, Ω p(x)dx = 0. 1
1877-0509 © 2017 The Authors. Published by Elsevier B.V. Peer-review under responsibility of the scientific committee of the International Conference on Computational Science 10.1016/j.procs.2017.05.093
1
2266
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
In the formulation (1), the velocity and pressure are treated as the two primal variables, so it is called two-field model. The three-field model of Stokes flow introduces the stress tensor as a new variable except for velocity and pressure. We define the stress tensor τxx τxy τ = . τyx τyy Then (1) can be reformulated as −∇ · τ + ∇p = f ,
η −1 τ = ∇u.
(3)
The three-field model has been extensively used in nonlinear Stokes problems [2, 6, 12]. The mixed finite element method is a very popular method for simulating flow problems in porous media [1,9], Stokes and Navier-Stokes equations [3,4,7]. In particular, a series of research activities is carried out on simple and efficient nonconforming finite element methods for the two-field model of Stokes problems, for example, [5, 10, 11]. These mixed methods discretize the velocity and pressure in different finite element spaces. On the other hand, the dual-mixed finite element methods for a three-field model of Stokes problems have also been developed, for example, [2,6,12]. Compared to mixed methods based on two-field model, dual-mixed methods for a three-field model provide more flexibility in approximating three variables using different finite element spaces, so they may achieve higher computational accuracy for not only the pair of velocity and pressure but also the stress. The relationship between mixed finite element methods and finite difference or finite volume schemes is an interesting research issue. It has been shown in [1] that the cell-centered finite difference method for elliptic problems can be derived from the lowest-order Raviart-Thomas (RT0 ) spaces [13] on rectangles using appropriate quadrature rules. The staggered grid introduced by Harlow and Welch [8] has been used extensively and popularly in computational fluid mechanics. Finite volume schemes on staggered grid have also been the basis of many industrial applications. It is shown in [7] that the MAC method can be derived by combining a velocity-vorticity mixed finite element method with an adequate quadrature formula. In this paper, we will first construct a three-field weak formulation for Stokes problems, and from this, a new dual-mixed finite element method is proposed. We discretize the velocity in the lowest-order Raviart-Thomas space and the pressure with piecewise constant functions, respectively. The components of stress tensor are approximated by piecewise constant functions or Q1 functions. Applying trapezoid quadrature rules to our proposed scheme, we derive a finite volume method on staggered grid, which is extensively used in computational fluid mechanics and engineering.
2
Finite volume method on staggered grid
We now formulate a finite volume method on staggered grid, which is one of the most common methods for simulating Stokes equations in practical applications. For sake of simplification, we focus on the case of rectangular domain, but the proposed scheme can be extended to the other regular domains. Moreover, the viscosity is assumed to be a constant. We consider a rectangular domain Ω = [0, Lx ] × [0, Ly ]. For sake of simplification, we use a uniform rectangular mesh for the domain, and denote the mesh size by h. The interval [0, Lx ] is uniformly divided into M subintervals, while the interval [0, Ly ] is uniformly divided into N subintervals. The variables located on the vertexes of each element are denoted by integer subscripts or superscripts, while variables on the middle points of an edge or element are denoted 2
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
by fraction subscripts or superscripts. For example, for a control volume [xi , xi+1 ] × [yj , yj+1 ], the pressure is located at its center (denoted by pi+ 12 ,j+ 12 ), the velocity component u is located at the middle of the left and right edges (denoted by ui,j+ 12 and ui+1,j+ 12 ), and the velocity components v is located at the middle of the top and bottom edges (denoted by vi+ 21 ,j and 1 1 1 1 i+ 1 ,j+ 1 T i+ 1 ,j+ 1 vi+ 21 ,j+1 ). We define the discrete source term f i+ 2 ,j+ 2 = fx 2 2 , fy 2 2 as f i+ 2 ,j+ 2 = xi+1 yj+1 1 f (x, y)dydx. The discrete form of Stokes equation on the interior elements is h2 x i yj expressed as η 4ui,j+ 12 − ui−1,j+ 12 − ui+1,j+ 21 − ui,j− 21 − ui,j+ 32 1 i− 1 ,j+ 1 i+ 1 ,j+ 1 (4) + h pi+ 21 ,j+ 12 − pi− 12 ,j+ 12 = h2 fx 2 2 + fx 2 2 , 2 η 4vi+ 21 ,j − vi− 12 ,j − vi+ 32 ,j − vi+ 21 ,j−1 − vi+ 21 ,j+1 1 i+ 1 ,j+ 1 i+ 1 ,j− 1 + h pi+ 12 ,j+ 12 − pi+ 12 ,j− 12 = h2 fy 2 2 + fy 2 2 . 2
(5)
The mass balance equation on the interior elements is discretized as ui+1,j+ 12 − ui,j+ 21 + vi+ 12 ,j+1 − vi+ 21 ,j = 0.
(6)
The discrete forms on the boundary elements can be similarly obtained by substituting the boundary conditions.
3
Dual-mixed finite elements for Stokes equation
We define a function space as 2 2×2 ∂ τˆxy ∂ τˆyx τˆxx τˆxy 2 2 ∈ L (Ω), ∈ L (Ω) . ∈ L (Ω) : S = τˆ = τˆyx τˆyy ∂x ∂y A bilinear operator is defined as ∂τyx ∂v ∂τxy ∂u τxx − v −u + τyy dx. A∇ (u, τ ) = ∂x ∂x ∂y ∂y Ω
(7)
A mixed weak form of Stokes equation is to seek (τ , p, u) ∈ S × L2 (Ω) × H(div; Ω) such that A∇ (v, τ ) − (p, ∇ · v) = (f , v), v ∈ H(div; Ω),
(8)
η −1 τ , τˆ = A∇ (u, τˆ ), τˆ ∈ S,
(9)
(∇ · u, w) = 0, w ∈ L2 (Ω).
(10)
We use a uniform rectangular mesh, which is the same to that stated in the above section. Let K be any element. The finite element space S h is defined as τ : τˆxx ∈ P0 (K) , τˆyy ∈ P0 (K) , τˆxy ∈ Q1 (K) , τˆyx ∈ Q1 (K)}. S h = {ˆ
(11) 3
2267
2268
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
The finite element space RT0 on an element K of a rectangular mesh is defined as � � RT0 = {v : v = a1K + a2K x, a3K + a4K y , aiK ∈ R, 1 ≤ i ≤ 4}.
(12)
More precisely, the basis function of RT0 can be expressed as v = [φxi,j+ 1 (x, y), 0]T or v = 2
[0, φyi+ 1 ,j (x, y)]T , where 2
φxi,j+ 1 (x, y) = 2
xi+1 −x xi+1 −xi , x−xi−1 xi −xi−1 ,
yj+1 −y yj+1 −yj , y−yj−1 yj −yj−1 ,
φyi+ 1 ,j (x, y) 2
=
0,
0,
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi−1 , xi ], y ∈ [yj , yj+1 ] , elsewhere
(13)
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi , xi+1 ], y ∈ [yj−1 , yj ] . elsewhere
(14)
Thus, the velocity uh ∈ RT0 has the unique representation uh (x) =
N M+1 ��
ui,j+ 12 φxi,j+ 1 (x), vh (x) = 2
i=1 j=1
+1 M N � �
vi+ 12 ,j φyi+ 1 ,j (x), x ∈ Ω.
The basis functions of the finite element space S h have the following forms � 1, x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] i+ 1 ,j+ 1 i+ 1 ,j+ 1 , τˆxx 2 2 (x, y) = τˆyy 2 2 (x, y) = 0, elsewhere
i,j i,j (x, y) = τˆxy (x, y) = τˆyx
xi+1 −x xi+1 −xi x−xi−1 xi −xi−1 xi+1 −x xi+1 −xi x−xi−1 xi −xi−1
yj+1 −y yj+1 −yj , yj+1 −y yj+1 −yj , y−yj−1 yj −yj−1 , y−yj−1 yj −yj−1 ,
0,
(15)
2
i=1 j=1
(16)
x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] x ∈ [xi−1 , xi ], y ∈ [yj , yj+1 ] x ∈ [xi , xi+1 ], y ∈ [yj−1 , yj ] , x ∈ [xi−1 , xi ], y ∈ [yj−1 , yj ] elsewhere
(17)
We define the finite element approximation of the stress as � � τxx,h τxy,h , x ∈ Ω, τ h (x) = τyx,h τyy,h
(18)
where τxx,h (x) =
N M � �
i+ 1 ,j+ 12
τxx 2
i+ 1 ,j+ 12
τˆxx 2
(x),
τyy,h (x) =
i=1 j=1
τxy,h (x) =
N M � �
i+ 1 ,j+ 12
τˆyy 2
(x),
(19)
i=1 j=1
+1 M+1 � N� i=1 j=1
i,j i,j τxy τˆxy (x),
τyx,h (x) =
+1 M+1 �N � i=1 j=1
We define the finite element space Wh on an element K as Wh = {w : w is constant in K}. 4
i+ 1 ,j+ 12
τyy 2
i,j i,j τyx τˆyx (x).
(20)
Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for Jisheng the three-field model
Kou and Sun
Furthermore, we define the basis functions of Wh as 1, x ∈ [xi , xi+1 ], y ∈ [yj , yj+1 ] . wi+ 21 ,j+ 12 (x) = 0, elsewhere
(21)
The discrete forms of the pressure and source term are expressed as ph (x) =
N M
pi+ 21 ,j+ 21 wi+ 12 ,j+ 21 (x), f h (x) =
N M
1
1
f i+ 2 ,j+ 2 wi+ 12 ,j+ 12 (x), x ∈ Ω. (22)
i=1 j=1
i=1 j=1
The mixed finite element formulation is to seek (τ h , ph , uh ) ∈ S h × Wh × RT0 such that A∇ (v, τ h ) − (ph , ∇ · v) = (f h , v), v ∈ RT0 ,
(23)
−1 η τ h , τˆ = A∇ (uh , τˆ ), τˆ ∈ S h ,
(24)
(∇ · uh , w) = 0, w ∈ Wh .
(25)
From the properties of mixed finite elements, the proposed scheme has approximate errors of at least one order, but because of limitation in the length of this article, the detailed analysis will be presented in our future work.
4
Equivalence forms
We now demonstrate the equivalence between the proposed finite element method with appropriate quadrature rules and the finite volume method on staggered grid stated in Section 2. For simplification, we derive this equivalence for the interior elements only. i+ 21 ,j+ 21 −1 0 τ ˆ xx , we Firstly, we consider the representation of η τ h , τˆ . Taking τˆ = 0 0 have −1 i+ 1 ,j+ 1 η τ h , τˆ = h2 η −1 τxx 2 2 . (26) 0 0 leads to The choice of τˆ = i+ 1 ,j+ 1 0 τˆyy 2 2
Let τˆ =
i,j 0 τˆxy 0 0
η −1 τ h , τˆ
−1 i+ 1 ,j+ 1 η τ h , τˆ = h2 η −1 τyy 2 2 .
(27)
, we obtain
−1 i,j η τxy,h , τˆxy xi+1 yj+1 i,j = η −1 τxy,h τˆxy dydx + =
xi
+
yj
xi+1
xi
yj
i,j η −1 τxy,h τˆxy dydx + yj−1
xi
xi−1 xi
xi−1
yj+1 i,j η −1 τxy,h τˆxy dydx
yj
yj i,j η −1 τxy,h τˆxy dydx.
(28)
yj−1
5
2269
2270
Dual-mixed finite elements for Jisheng the three-field model Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274
Kou and Sun
We apply trapezoid quadrature rule to the first term on the right-hand side of (28) � xi+1 � yj+1 i,j η −1 τxy,h τˆxy dydx xi
=
�
yj
xi+1 �
xi
=
≈
yj+1
yj
j+1 i+1 � � k,l k,l i,j η −1 τxy τˆxy τˆxy dydx k=i l=j
�2 yj+1 − y dx dy yj+1 − yj xi yj �2 � xi+1 � yj+1 � xi+1 − x x − xi yj+1 − y −1 i+1,j dx dy +η τxy xi+1 − xi xi+1 − xi yj+1 − yj xi yj �2 � yj+1 � xi+1 � yj+1 − y y − yj xi+1 − x i,j+1 dx dy +η −1 τxy x − x y i+1 i j+1 − yj yj+1 − yj yj xi � xi+1 � yj+1 yj+1 − y y − yj xi+1 − x x − xi −1 i+1,j+1 dx dy +η τxy x − x x − x y i+1 i i+1 i j+1 − yj yj+1 − yj yj xi i,j η −1 τxy
�
xi+1
�
xi+1 − x xi+1 − xi
�2
�
yj+1
�
h2 i,j τ . 4η xy
(29)
The rest three terms on the right-hand side of (28) can be approximated by the similar way, and consequently, in this case, we get � −1 � � � i,j i,j η τ h , τˆ ≈ h2 η −1 τxy = η −1 τxy,h , τˆxy . (30) By the similar appraoches, we take τˆ =
�
0 i,j τˆyx
0 0
�
and obtian
� � � i,j i,j η −1 τ h , τˆ = η −1 τyx,h , τˆyx ≈ h2 η −1 τyx .
�
Secondly, we consider the form of A∇ (uh , τˆ ). We take � � i+ 1 ,j+ 1 � τˆxx τˆxy τˆxx 2 2 = τˆ = τˆyx τˆyy 0
0 0
�
(31)
,
and then A∇ (uh , τˆ ) becomes A∇ (uh , τˆ )
Similarly, let τˆ =
�
� � � xi+1 � yj+1 � ∂uh ∂uh i+ 12 ,j+ 21 τˆxx dx = τˆxx dx ∂x ∂x Ω xi yj � xi+1 � yj+1 � � ∂uh dx = h ui+1,j+ 21 − ui,j+ 12 . = ∂x xi yj
=
0
� �
�
0 i+ 21 ,j+ 12
0 τˆyy
A∇ (uh , τˆ ) =
�
6
, we get
xi+1
xi
(32)
�
yj+1
yj
� � ∂vh dx = h vi+ 12 ,j+1 − vi+ 12 ,j . ∂x
(33)
Dual-mixed finite elements for Jisheng the three-field model Kou et al.Stokes / Procedia Computer Science 108C (2017) 2265–2274
Taking τˆ =
i,j 0 τˆxy 0 0
Kou and Sun
gives
i,j ∂ τˆxy dx −vh ∂x Ω xi+1 yj+1 xi yj+1 i,j i,j ∂ τˆxy ∂ τˆxy − dydx − dydx vh vh ∂x ∂x xi yj xi−1 yj xi yj xi+1 yj i,j i,j ∂ τˆxy ∂ τˆxy dydx − dydx − vh vh ∂x ∂x xi yj−1 xi−1 yj−1 1 yj+1 yj+1 − y y − yj (yj+1 − y) dy vi+ 21 ,j + vi+ 21 ,j+1 h yj yj+1 − yj yj+1 − yj 1 yj+1 yj+1 − y y − yj vi− 21 ,j (yj+1 − y) dy − + vi− 21 ,j+1 h yj yj+1 − yj yj+1 − yj 1 yj yj − y y − yj−1 vi+ 21 ,j−1 (y − yj−1 ) dy + + vi+ 21 ,j h yj−1 yj − yj−1 yj − yj−1 1 yj yj − y y − yj−1 (y − yj−1 ) dy. v 1 − + vi− 21 ,j h yj−1 i− 2 ,j−1 yj − yj−1 yj − yj−1
A∇ (uh , τˆ ) = =
=
(34)
Applying trapezoid quadrature rule to the terms on the right-hand side of (34), we deduce A∇ (uh , τˆ ) ≈ h vi+ 21 ,j − vi− 12 ,j . (35) Taking τˆ = deduce
0 i,j τˆyx
0 0
, and furthermore, by the similar approaches used in (34)-(35), we A∇ (uh , τˆ ) ≈ h ui,j+ 21 − ui,j− 12 .
(36)
Combining the above derivations, we derive the relations between τ h and uh as η η i+ 1 ,j+ 1 i+ 1 ,j+ 1 ui+1,j+ 12 − ui,j+ 12 , τyy 2 2 = vi+ 12 ,j+1 − vi+ 21 ,j , τxx 2 2 = h h η ui,j+ 12 − ui,j− 21 . h T Thirdly, we consider the expressions of A∇ (v, τ h ). Take v = φxi,j+ 1 , 0 i,j ≈ τxy
η vi+ 12 ,j − vi− 21 ,j , h
i,j τyx ≈
(37)
(38)
2
A∇ (v, τ h )
=
∂φx 1 i,j+ 2
Ω
∂x
τxx,h −
φxi,j+ 1 2
∂τyx,h ∂y
dx.
(39)
Using their definitions, the first term on the right-hand side of (39) is calculated as ∂φx 1 xi+1 yj+1 ∂φx 1 i,j+ 2 i+ 12 ,j+ 12 i,j+ 2 τxx,h dx = τxx dydx ∂x ∂x xi yj Ω 7
2271
2272
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model
+
xi
xi−1
=
yj+1
yj
∂φxi,j+ 1
i− 1 ,j+ 1 τxx 2 2
2
∂x
Kou and Sun
dydx
i+ 1 ,j+ 1 i− 1 ,j+ 1 −h τxx 2 2 − τxx 2 2 .
(40)
The second term on the right-hand side of (39) is rewritten as xi+1 yj+1 xi yj+1 ∂τyx,h ∂τyx,h ∂τyx,h dx = dydx + dydx φxi,j+ 1 φxi,j+ 1 φxi,j+ 1 2 2 2 ∂y ∂y ∂y Ω xi yj xi−1 yj xi+1 i,j+1 xi+1 − x xi+1 − x i,j = dx τyx − τyx xi+1 − xi xi+1 − xi xi xi+1 i+1,j+1 xi+1 − x x − xi i+1,j dx − τyx τyx + xi+1 − xi xi+1 − xi x xi i x − xi−1 i,j+1 x − xi−1 i,j dx τyx − τyx + x − x x i i−1 i − xi−1 xi−1 xi i−1,j+1 x − xi−1 xi − x i−1,j + dx. (41) τyx − τyx xi − xi−1 xi − xi−1 xi−1 Applying trapezoid quadrature rule to the four terms on the right-hand side of (41) leads to i,j+1 ∂τyx,h i,j dx ≈ −h τyx −φxi,j+ 1 . (42) − τyx 2 ∂y Ω T Let v = 0, φyi+ 1 ,j , we have 2
A∇ (v, τ h ) =
∂φy 1 i+ ,j 2
∂y
Ω
τyy,h −
φyi+ 1 ,j 2
∂τxy,h ∂x
(43)
dx.
The first term on the right-hand side of (43) is directly calculated as ∂φy 1 i+ ,j 2
Ω
∂y
τyy,h dx =
xi+1
yj
xi+1
xi
+
2
yj−1
xi
=
∂φyi+ 1 ,j yj+1
yj
∂y
i+ 1 ,j− 12
τyy 2
∂φyi+ 1 ,j 2
∂y
dydx
i+ 1 ,j+ 12
τyy 2
i+ 1 ,j+ 1 i+ 1 ,j− 1 −h τyy 2 2 − τyy 2 2 .
dydx (44)
For the second term on the right-hand side of (43), we use the trapezoid quadrature rule xi+1 yj xi+1 yj+1 ∂τxy,h ∂τxy,h ∂τxy,h y y − dx = − dydx − dydx φi+ 1 ,j φyi+ 1 ,j φi+ 1 ,j 2 2 2 ∂x ∂x ∂x xi yj−1 xi yj Ω yj+1 i+1,j yj+1 − y yj+1 − y i,j = − dy τxy − τxy yj+1 − yj yj+1 − yj yj yj+1 i+1,j+1 yj+1 − y y − yj i,j+1 − dy τxy − τxy yj+1 − yj yj+1 − yj yj 8
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model yj
Kou and Sun
y − yj−1 y − yj−1 dy yj − yj−1 yj − yj−1 yj−1 yj i+1,j−1 y − yj−1 yj − y i,j−1 − dy τxy − τxy yj − yj−1 yj − yj−1 yj−1 i+1,j i,j ≈ −h τxy − τxy . (45) T Finally, we consider the formulations of (ph , ∇ · v) and (f h , v). Taking v = φxi,j+ 1 , 0 2 T y and v = 0, φi+ 1 ,j respectively, we obtain 2 ∂φxi,j+ 1 2 = −h pi+ 12 ,j+ 12 − pi− 21 ,j+ 12 , (46) (ph , ∇ · v) = ph , ∂x −
(ph , ∇ · v) =
ph ,
i,j i+1,j − τxy τxy
∂φyi+ 1 ,j 2
∂y
= −h pi+ 12 ,j+ 12 − pi+ 12 ,j− 12 .
(47)
T T respectively, we get For (f h , v), taking v = φxi,j+ 1 , 0 and v = 0, φyi+ 1 ,j 2
2
1 i+ 1 ,j+ 1 i− 1 ,j+ 1 (f h , v) = fx,h , φxi,j+ 1 = h2 fx 2 2 + fx 2 2 , 2 2
(48)
1 i+ 1 ,j+ 1 i+ 1 ,j− 1 (f h , v) = fy,h φyi+ 1 ,j = h2 fy 2 2 + fy 2 2 . 2 2
(49)
Thus the equation (24) can be rewritten as i− 1 ,j+ 1 i+ 1 ,j+ 1 i,j+1 i,j −h τxx 2 2 − τxx 2 2 + τyx + h pi+ 21 ,j+ 12 − pi− 12 ,j+ 21 − τyx 1 i+ 1 ,j+ 1 i− 1 ,j+ 1 = h2 f x 2 2 + f x 2 2 , 2 i+ 1 ,j+ 1 i+ 1 ,j− 1 i+1,j i,j −h τyy 2 2 − τyy 2 2 + τxy + h pi+ 21 ,j+ 12 − pi+ 12 ,j− 12 − τxy 1 i+ 1 ,j+ 1 i+ 1 ,j− 1 = h2 f y 2 2 + f y 2 2 . 2
(50)
(51)
The discrete forms (4) and (5) are obtained substituting (37) and (38) into the equations (50) and (51). For the equation (∇ · uh , w) = 0, let w = 1 on an element K, we get (6). We note that the above derivations are carried out on the interior elements, but the boundary elements can be handled by the similar approaches and taking into account the boundary conditions.
5
Numerical results
We test the benchmark driven cavity problem. For Ω = [0, 1]2 , we have that uB = [1, 0]T on the top boundary (that is, the portion of the boundary satisfying 0 ≤ x ≤ 1 and y = 1) and uB = 0 on the rest portions of the boundary. The source term is zero. The viscosity is set as η = 1. Figs. 1 shows the computed flow quiver on a rectangular mesh composed of 30 × 30 square elements. 9
2273
Jisheng Kou et al. / Procedia Computer Science 108C (2017) 2265–2274 Dual-mixed finite elements for the three-field Stokes model
Kou and Sun
Flow quiver
1 0.8 0.6
Y
2274
0.4 0.2 0 0
0.2
0.4
X
0.6
0.8
1
Figure 1: The flow quiver.
6
Conclusions
A new weak formulation is developed for the three-field model of Stokes problems. Based on this formulation, we proposed a dual-mixed finite element method on a rectangular mesh. A finite volume method on staggered grid is derived by applying trapezoid quadrature rules for the proposed scheme.
References [1] T. Arbogast, M.F. Wheeler and I. Yotov. Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences. SIAM J. Numer. Anal., 34: 828–852, 1997. [2] J. Baranger, K. Najib and D. Sandri. Numerical analysis of a three-beads model for a quasiNewtonian flow. Comput. Methods Appl. Mech. Engrg., 109: 281–292, 1993. [3] S. C. Brenner and R. Scott. The mathematical theory of finite element methods, volume 15. Springer, 2007. [4] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer Verlag, 1991. [5] M. Crouzeix and P.A. Raviart. Conforming and nonconforming finite element methods for solving the stationary stokes equations. RAIRO Analysis Numerical, 7:33–36, 1973. [6] V. J. Ervin, J. S. Howell and I. Stanculescu. A dual-mixed approximation method for a threefield model of a nonlinear generalized Stokes problem. Comput. Methods Appl. Mech. Engrg., 197: 2886–2900, 2008. [7] V. Girault, H. Lopez. Finite-element error estimates for the MAC scheme. IMA Journal of Numerical Analysis, 16(3): 347-379, 1996. [8] F. H. Harlow and J. E. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys. Fluids, 8: 2182–2189, 1965. [9] J. Kou, S. Sun and Y. Wu. Mixed finite element-based fully conservative methods for simulating wormhole propagation. Comput. Methods Appl. Mech. Engrg., 298: 279–302, 2016. [10] B. P. Lamichhane. A nonconforming finite element method for the Stokes equations using the Crouzeix-Raviart element for the velocity and the standard linear element for the pressure. Int. J. Numer. Meth. Fluids, 74: 222–228, 2014. [11] J. Li and Z. Chen. A new local stabilized nonconforming finite element method for the Stokes equations. Computing, 82:157–170, 2008. [12] H. Manouzi and M. Farhloul. Mixed finite element analysis of a nonlinear three-fields Stokes model. IMA J. Numer. Anal., 21: 143–164, 2001. [13] R.A. Raviart and J.M. Thomas. Primal hybrid finite element methods for 2nd order elliptic equations. Math. Comput., 31: 391–413, 1977.
10