# A Two-Level Nonconforming Rotated Quadrilateral Finite Element Method for the Stationary Navier–Stokes Equations

Authors: Weijun Tian; Liquan Mei; Yinnian He
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we propose a two-level nonconforming rotated finite element (TNRFE) method for solving the Navier–Stokes equations. A new nonconforming rotated finite element (NRFE) method was proposed by Douglas added by conforming bubbles to velocity and discontinuous piecewise constant to the pressure on quadrilateral elements possessing favorable stability properties. The TNRFE method involves solving a small Navier–Stokes problem on a coarse mesh with mesh sizeH and a large linearized Navier–Stokes problem on a fine mesh with mesh size h by the NRFE method. If we choose h=OH2, the TNRFE method gives the convergence rate of the same order as that of the NRFE method. Compared with the NRFE method, the TNRFE method can save a large amount of CPU time. In this paper, the stability of the approximate solutions and the error estimates are proved. Finally, the numerical experiments are given, and results indicate that the method is practicable and effective.


---

## Body

## 1. Introduction

When the Navier–Stokes equations are discretized by finite element methods, two problems often arise: one is that the discrete inf-sup condition is broken, and the other is that the pseudo-oscillation is caused by the dominant convection term. In order to solve the first problem, we use the new nonconforming element proposed by Jim Douglas [1] added by conforming bubbles to the velocity and discontinuous piecewise constant to the pressure on quadrilateral elements which possess favorable stability properties. Brezzi and Russo [2] found that adding and eliminating bubbles to the finite element space is equivalent to the addition of a stabilizing term of a streamline diffusion type. It is equivalent to tuning up the amount of streamline artificial viscosity and corresponds to reducing the residual inside each element.To deal with the second problem, we adopt the Newton TNRFE method which involves solving one small Navier–Stokes problem on a coarse mesh with mesh sizeH and a large linearized Navier–Stokes problem on a fine mesh with mesh size h.It is well known that the approximate space of velocity and pressure should generally satisfy the LBB condition when solving Navier–Stokes equations by finite element methods. This matching is difficult for the conforming element. For example, if the triangle element is used, the velocity is a piecewise linear element, and the pressure is a piecewise constant, then the LBB condition is not satisfied. If the velocity is changed to piecewise quadratic, then the error estimates for velocity will lose one order. If you use the following two rectangular elements, (1) the velocity is a bilinear polynomial interpolation and the pressure is a piecewise constant; (2) the velocity is double polynomial interpolation and the pressure is of bilinear form, then the LBB condition is not satisfied.To some extent, the nonconforming element can overcome this difficulty, and it is often used in practical computation because of its simple structure, economical calculation, and error matching. The nonconforming finite element method for the Stokes problem was first proposed by Crouzeix and Raviart [3], and they used the piecewise linear triangle element with three midpoints as the velocity approximation space and the piecewise constant finite element as the pressure approximation space to obtain an approximation scheme (the C-R scheme) for the Stokes problem. In this way, not only the LBB condition is satisfied but also some optimal error estimates of velocity and pressure are obtained. Temam [4] also proposed the nonconforming linear element as the approximate space of velocity, which obtained satisfactory error estimates. Because the rectangular element has the advantages of simple mesh and less bandwidth of the stiffness matrix, Rannacher and Turek [5] analyzed the rotational bilinear nonconforming element with local base span 1,x,y,x2−y2. However, if the definition of the global nonconforming space which requires continuity at the midpoint of the common inner boundary of the adjacent elements is adopted, the optimal error estimates for the real quadrilateral subdivision region will not be obtained. To solve this problem, Jim Douglas et al. [1] modified the rotated bilinear local basis span 1,x,y,x2−y2 to span 1,x,y,3x2−5x4−3y2−5y4, which has the following properties: <1,wj−wk>Γjk=0 and <1,wj>Γj=0. Han [6] proposed that a nonconforming rectangular element with five degrees of freedom for the velocity, whose node parameters are the function value at the center and the midpoint of the four sides of the element, and the shape function space span 1,x,y,1/25x4−3x2,1/25y4−3y2 for the standard element −1,12 can also obtain the best order error estimate. Although the element has the same order of convergence speed as the Douglas rotation element, it has one more degree of freedom than the latter, which increases the computational load and decreases the efficiency. In solving the Stokes problem and the Navier–Stokes problem, Cai et al. [7, 8] used the Douglas rotation element, added a compatible bubble function to the velocity approximation space as the internal degree of freedom, and used the piecewise constant finite element space as the pressure approximation space. Cai et al. [7] gave the existence and uniqueness of the solution and the convergence rate of the nonconforming Galerkin (NFG) method when the Douglas rotation element was used for the velocity, but the numerical experiments were not given.A part of the work of this paper is to solve the steady Stokes equations and Navier–Stokes equations by using the NRFE method, explain the existence and uniqueness of the solution and the convergence speed, and give the concrete design process and numerical experiment of the algorithm.In practice, solving the Navier–Stokes equations is very time-consuming, so people try to save computing time as much as possible. One of them is the two-layer grid method, which is to solve a nonlinear equation on a coarse grid, and the purpose of solving a linear equation on a fine mesh is to obtain the approximate solution of the nonlinear equation in less time and to maintain the optimal convergence speed, specifically. The two-layer mesh method is based on the finite element spaceXH,Xh for the velocity approximation (H is the coarse mesh size parameter, and h is the fine mesh size parameter, h≪H) and on the finite element space MH,Mh for the pressure approximation; the first step is to get the solution of uH,pH on a coarse grid by the NRFE method. Generally speaking, because the grid is coarse, it does not take much time. The second step is to solve a linear equation using known uH,pH. In some papers, the third step is to solve a correction problem on the coarse grid by the correction technique. In recent years, the two-layer grid method has attracted more and more attention. Layton [9], Li and Hou [10], He and Li [11], Mei [13], Shi and Wang [14], Chen [15], and others have done a lot of work in this field.The other main work of this paper is to prove the stability and convergence rate of the TNRFE method for solving Navier–Stokes equations. The design process and numerical experiment of the algorithm are given.The TNRFE method is used to solve a nonlinear problem on a coarse grid with a grid size ofH and a linear problem on a fine grid with a grid size of h. However, the NRFE method is used to solve a nonlinear problem on a fine grid with a grid size of h. Therefore, the TNRFE method is simpler than the NRFE method in calculation and can save a lot of CPU time, and the TNRFE method has the same order of convergence speed as the NRFE method when the grid size is selected properly. In this paper, the numerical experiments and theoretical analysis prove that the Newton TNRFE method and the NRFE method have the same order of convergence speed under the condition of h=OH2.In Section2, we briefly review some knowledge and important properties of Sobolev spaces and give the variational forms and the existence, uniqueness, and regularity estimates of solutions for the Navier–Stokes equations. In Section 3, we first introduce the construction and properties of the nonconforming rectangular element space. We explain the NRFE method, existence, uniqueness, and convergence rate of the approximation solution for Stokes equations. The design process and finite element analysis of the algorithm are given in detail. In Section 4, the NRFE method, existence, uniqueness, and convergence rate of the approximate solution for Navier–Stokes equations are discussed. In Section 5, the stability and convergence rate of the solution of the TNRFE method for Navier–Stokes equations are proved. In Section 6, a large number of numerical experiments are given to verify the correctness and effectiveness of the algorithm.


## 2. Governing Equations

In this paper, we study the linearized energy-conservative finite element method for the following nonlinear Navier–Stokes equations:(1)−νΔu+u⋅∇u+∇p=f,∀x∈Ω,∇⋅u=0,∀x∈Ω,u=0,∀x∈Γ,where Ω is a bounded domain and convex domain in Rdd=2,3, u=u1x,u2x is the velocity vector, p=px is the pressure, f=fx is the external force density, and ν>0 is the dynamic viscosity coefficient.In order to describe the variational form of Navier–Stokes equations, the Sobolev space is introduced:(2)X=H01Ω2,Y=L2Ω2,M=L02Ω=q∈L2Ω;∫Ωqdx=0,where ⋅,⋅,‖⋅‖0 denotes the inner product and norm on L2Ωd and X is equipped with the usual inner product and norm.(3)u,v=∇u,∇v,∇u0=u,u1/2,and ⋅,⋅Γ,‖|⋅‖|Γ denotes the inner product and norm over L2Γd.Define, respectively, the continuous bilinear form onX×X and X×M:(4)au,v=νu,v,∀u,v∈X,dv,q=∇⋅v,q,∀v∈X,q∈M.In addition, we introduce a closed subset ofX:(5)V=v∈X;dv,q=0,∀q∈M.LetAu=−Δu. As we all know, A is a linear, unbounded, and self-conjugate operator in Y, defined in(6)DA=H2Ω2∩X.Theorem 1 (see [4]).
LetΩ be a smooth or convex polygon of C2; for given f∈Y, the Stokes problem(7)−νΔu+∇p=f,∀x∈Ω,∇⋅u=0,∀x∈Ω,u=0,∀x∈Γ,there exists a unique solution of u,p∈DA×H1Ω∩M and satisfies(8)‖u‖2+‖p‖1≤c‖f‖0,where c is a positive constant dependent on Ω and ‖⋅‖i is a norm over HiΩ or HiΩ2i=1,2.Define the trilinear form onX×X×X:(9)bu,v,w=u⋅∇v,w+12∇⋅uv,w=12u⋅∇v,w−12u⋅∇w,v,∀u,v,w∈X,which satisfies the following properties:(10)u⋅∇v,w+u⋅∇w,v=∫Γw⋅vu⋅nds−∫Ωw⋅v∇⋅udx,∀u,v,w∈X|bu,v,w|≤Cb‖∇u‖0‖∇v‖0‖∇w‖0,∀u,v,w∈Xbu,v,w=−bu,w,v,∀u,v,w∈X.So, the variational form of problem (1) is(11)findu,p∈X×M,soau,v−dv,p+bu,u,v=f,v,∀v∈X,du,q=0,∀q∈M,and problem (2) in the variational form is(12)findu,p∈X×M,so,au,v−dv,p=f,v,∀v∈X,du,q=0,∀q∈M.The inequalities that will be used in this article are described here first.(1)
Hölder’s inequality:Ifu∈LpΩandv∈LqΩ,p,q are positive real numbers, then(13)∫ΩuvdxuLpvLq,1p+1q=1.(2)
Young’s inequality:Ifa,b,p,q,andϵ are all positive real numbers, then(14)ab≤ϵapp+ϵ−q/pbqq,1p+1q=1.(3)
Sobolev’s embedding theorem:Letm>0 be an integer, 1≤p<∞, and let Ω be an open subset of Rn and have Lipschitz continuous boundary, so the following embedding relationship holds:(15)Wm,pΩ⟶LqΩ,q=npn−mp,mp<nLqΩ,∀q∈1,∞,mp=nC0,m−n/pΩ¯,mp>n>m−1pC0,αΩ¯,∀α∈0,1,n=m−1pC0,1Ω¯,n<m−1p.In this paper, we give the existence and uniqueness of solutions for the classical Navier–Stokes equations.Theorem 2 (see [16]).
Ifv and f∈X′ satisfy(16)Cbν2f−1<1,f−1=supv∈Xf,v∇v0,then problem (11) has a unique solution of u,p∈DA×H1Ω∩M, and u∈V satisfies(17)∇u0≤1νf−1.Theorem 3 (see [12]).
Iff∈Y satisfies uniqueness condition (8), then the solution of problem (11) satisfies the following regularity estimate:(18)u2+p1≤C,where C is a positive constant dependent on ν,f,Ω.


## 3. A Nonconforming Rotated Finite Element Method for the Stokes Equations

In the first section of this section, we review the construction and properties of the nonconforming rotated rectangular element space. In the second section, we discuss the existence, uniqueness, and stability of the solution of the problem approximated by the NRFE method. In the third section, we give the convergence rate of this method. In the fourth section, the concrete design of the algorithm is introduced. The key of solving the Stokes equations with the NRFE method is to condense the internal degree of freedom in the process of establishing the stiffness. In the fifth section, we describe the finite element analysis of the method.In this paper, we use a rectangular partition of the same scale for the velocity and pressure. Here and later,h is a positive real parameter that tends to zero, and 0<h≤1.Definition 1.
Letτh=τhΩ be a partition of the region Ω¯ and Ωj be any unit of τh; if τh satisfies the following conditions:(i)
There exists a constantσ>0 so that hj/ρj≤σ,∀Ωj∈τh, where hj=diamΩj=sup|x−y|,∀x,y∈Ωjandρj=supdiamS,S is any inside ball of Ωj(ii)
There exists a constantγ>0 that has nothing to do with h such that maxh/hj,∀Ωj∈τh≤γ, where h=maxΩjhjthenτh is called the quasi-uniform regular partition of the region Ω¯.
LetΩ¯=∪j=1JΩ¯j⊂R2,andΩj is any quasi-uniform regular quadrilateral element of the region Ω¯, recorded as(19)Γj=Γ∩Γj,Γjk=Γkj=∂Ωj∩Γk.
The midpoints ofΓj and Γjk are, respectively, ξj and ξjk. We define the standard base on the standard unit:(20)QK^=Span1,ξ,η,ξη,ξ2−53ξ4−η2−53η4,(21)QΩj=v:v=v^°Fj−1,v^∈QK^,where Fj is a reversible affine transformation and FjK^=Ωj.
Now, we construct the finite element space.(1)
Nonconforming rectangular finite element space(22)Xh=v:vj=v|Ωj∈QΩj×QΩj,vjξjk=vkξkj,vjξj=0,∀j,ksatisfies the following properties:(i)
Inverse inequality:(23)‖∇χ‖0,j≤chj−1‖χ‖0,j,∀χ∈Xh.(ii)
Poincaré–Friedrichs inequality:(24)‖v‖0,j2≤h2‖∇v‖0,j2+∫Γ1vds2,∀v∈H1Ωj2,Γ1⊆∂Ωj.(iii)
Trace inequality:(25)‖|v‖|j≤chj−1/2‖v‖0,j+‖v‖0,j1/2‖∇v‖0,j1/2,∀v∈H1Ωj.2(iv)
(26)∫Γjkvids=0,∫Γjvids=0,∀v=v1,v2T∈Xh,i=1,2,where‖⋅‖m,j2=‖⋅‖HmΩj22, ⋅j2=∑k∈Zj⋅L2Γjk22,andZj=k,Ω¯k∩Ω¯j≠∅. Γjk can be replaced by Γj at the outer boundary, and vi=vi|Γjk−vi|Γkj means the jump of vi on the public boundary Γjk, 1≤j≤J.(2)
Piecewise constant finite element space(27)Mh=q∈M:q|Ωj∈P0Ωj,∀jsatisfies the following properties:(i)
Inverse inequality:(28)‖∇q‖0,j≤chj−1‖q‖0,j,∀q∈Mh.(ii)
Poincaré–Friedrichs inequality:(29)‖p‖0,j2≤h2‖∇p‖0,j2+∫Γ1pds2,∀p∈H1Ωj,Γ1⊆∂Ωj.(iii)
Trace inequality:(30)‖|p‖|j≤chj−1/2‖p‖0,j+‖p‖0,j1/2‖∇p‖0,j1/2,∀p∈H1Ωj,where‖⋅‖m,j2=‖⋅‖HmΩj2, ⋅j2=∑k∈Zj⋅L2Γjk2,andZj=k,Ω¯k∩Ω¯j≠∅. Γjk can be replaced by Γj at the outer boundary, 1≤j≤J.
Let⋅,⋅j=⋅,⋅Ωj,1≤j≤J, define the discrete bilinear form(31)ahu,v=ν∑j=1J∇u,∇vj,∀u,v∈X∪Xhdhv,q=∑j=1J∇⋅v,qj,∀v∈Xh,q∈Mh,and norm(32)‖v‖1,h=∑j=1J‖∇v‖0,j21/2=1νahv,v,∀v∈Xh.
In addition, we introduce a closed subset ofXh:(33)Vh=v∈Xh;dhv,q=0,∀q∈Mh.
Define operatorπj:H1Ωj⟶QΩj, which satisfies(34)∫ejiπjvds=∫ejivds,∀v∈H1Ωj,i=1,2,3,4,and is equivalent to(35)πjvaji=1eji∫ejivds,i=1,2,3,4,where ejii=1,2,3,4 are the four edges of Ωj, and their corresponding midpoints are ajii=1,2,3,4, satisfying the following properties:(i)
(36)‖∇v−πjv‖0,j≤c‖∇v‖0,j,∀v∈H1Ωj.(ii)
(37)‖∇πjv‖0,j≤c‖∇v‖0,j,∀v∈H1Ωj.
Define operatorΠh:X⟶Xh, which satisfies(38)Πhvi|Ωj=πjvi,∀v∈X,∀j,i=1,2,and satisfies the following properties:(i)
(39)dhΠhv−v,q=0,∀q∈Mh.(ii)
(40)‖Πhv‖1,h≤c‖∇v‖0,∀v∈X.
Define the projection operatorRh:H2Ω2⟶Xh, which satisfies(41)Rhvξ=vξ,∀ξ=ξjkorξj,and P0:H1Ω2⟶Λh×Λh satisfies(42)P0vj,zΓ=∂vj∂nj,zΓ,∀z∈P0Γ2,∀Γ=ΓjkorΓj,where if record trτ=τ:δ,δ=1001,andσ:τ=∑i=12∑j=12σijτij, then(43)Λh=λ:λjk=trΓjkλ|Ωj∈P0Γjk;λjk+λkj=0;λj=trΓjλ|Ωj∈P0Γj.
Define the projection operatorSh:H1Ω⟶Mh, which satisfies(44)Shq,z=q,z,∀z∈Mh,and Q0:H1Ω⟶P0Γ satisfies(45)Q0q,zΓ=q,zΓ,∀z∈P0Γ,∀Γ=ΓjkorΓj.
They satisfy the following properties:(i)
Orthogonality:(46)P0vj,wjΓjk+P0vk,wkΓkj=P0vj,wj−wkΓjk=0,∀w∈Xh.(ii)
Approximation property:(47)‖v−Rhv‖0+h∑j=1J‖∇v−Rhv‖0,j21/2+h1/2∑j=1J‖|v−Rhv‖|j21/2+h3/2∑j=1J‖|∂v∂nj−P0v‖|j21/2≤ch2‖v‖2,v∈H2Ω2‖q−Shq‖0+h1/2∑j=1J‖|q−Q0q‖|j21/2≤ch‖q‖1,q∈H1Ω.
Here and later,nj,1≤j≤J, represents the unit outer normal vector of ∂Ωj.
Then, the nonconforming finite element approximation problem of (12) is(48)finduh,ph∈Xh×Mh,soahuh,v−dhv,ph=f,v,∀v∈Xh,dhuh,q=0,∀q∈Mh.Proposition 1 (see [7]).
If the bilinear formd⋅,⋅ satisfies the inf−sup condition, there is a constant β>0 such that(49)supv∈Xdv,q‖∇v‖0≥β‖q‖0,∀q∈M,then the bilinear form dh⋅,⋅ satisfies the discrete inf−sup condition: there is a positive constant β0 such that(50)supv∈Xhdhv,q‖v‖1,h≥β0‖q‖0,∀q∈Mh.Theorem 4 (see [7]).
Under the condition of Proposition1, if Vh is not empty, problem (48) has a unique solution uh,ph∈Xh×Mh, and uh∈Vh satisfies(51)‖uh‖1,h≤1ν‖f‖∗,‖f‖∗=supv∈Xhf,v‖v‖1,h.Proposition 2 (see [7]).
Letu,p and uh,ph be, respectively, the solutions of problems (12) and (48); there is a positive constant of c such that(52)‖u−uh‖1,h≤cinfv∈Xh‖u−v‖1,h+supv∈Vhahu,v−f,v‖v‖1,h‖p−ph‖0≤cinfq∈Mh‖p−q‖0+supv∈Xhahu,v−dhv,p−f,v‖v‖1,h+infv∈Xh‖u−v‖1,h+supv∈Vhahu,v−f,v‖v‖1,h.Proposition 3 (see [7]).
Letf,gj=∫∂Ωjfgds, for ∀ϕ,w∈X∪Xh, satisfy the following properties:(53)|∑j=1J∂w∂nj,ϕj|≤ch‖w‖2‖ϕ‖1,h,∀w∈DA|∑j=1Jq,ϕ⋅njj|≤ch‖q‖1‖ϕ‖1,h,∀q∈H1Ω.Theorem 5 (see [7]).
Letu,p∈H2Ω2×H1Ω and uh,ph∈Xh×Mh be, respectively, the solutions of problems (12) and (48); then,(54)‖u−uh‖1,h+‖p−ph‖0≤ch‖u‖2+‖p‖1.Using the dual argument method, we estimate the error ofL2.Consider the linear duality problem:(55)−νΔψ+∇χ=u−uh∀x∈Ω,∇⋅ψ=0∀x∈Ω,ψ=0∀x∈Γ.The variational form is(56)findψ,χ∈X×M,soav,ψ−dv,χ=u−uh,v,∀v∈X,dψ,q=0,∀q∈M.Ifu is the nonsingular solution of problem (2), then problem (56) has a unique solution (see [16]).In order to establish an estimate of theL2 error of velocity, Aubin–Nitsche argumentation method is used here. Let problem (34) be H2 regular, i.e.,(57)‖ψ‖2+‖χ‖1≤c‖u−uh‖0.Letψ,χ be the solution of problem (56), and ψh,χh∈Xh×Mh, satisfying(58)‖ψ−ψh‖1,h+‖χ−χh‖0≤ch‖ψ‖2+‖χ‖1.Theorem 6 (see [7]).
Letu,p∈H2Ω2×H1Ω and uh,ph∈Xh×Mh be, respectively, the solutions of problems (12) and (48); if problem (55) is H2 regular, then(59)‖u−uh‖0≤ch2‖u‖2+‖p‖1.
This section presents the algorithm design for solving problem (48).


### 3.1. Step 1: Condensation of Internal Degrees of Freedom

Let the four edges of the standard cellK^=−1,12 be, respectively, e^1,e^2,e^3,ande^4; the corresponding midpoints are a^1=0,−1,a^2=1,0,a^3=0,1,anda^4=−1,0, and the shape functions for each node of the Douglas rotation element obtained are, respectively,(60)φ1ξ,η=14−12η+38ξ2−53ξ4−η2−53η4φ2ξ,η=14+12ξ−38ξ2−53ξ4−η2−53η4φ3ξ,η=14+12η+38ξ2−53ξ4−η2−53η4φ4ξ,η=14−12ξ−38ξ2−53ξ4−η2−53η4.It is easy to verify(61)∫e^iφjds^=δij|e^i|,i,j=1,2,3,4,where δij is the Kronecker symbol and e^i is the length of the edge of e^i.The element speed interpolation with additional node-free speed terms can be expressed as follows:(62)u1=∑i=14φiui1+λ1ξη,u2=∑i=14φiui2+λ2ξη,where λ1andλ2 are the undetermined coefficient in the element’s internal definition, known as the internal degree of freedom.The expression is expressed in the matrix form:(63)u=ua+uλ,i.e.,(64)u1u2=φ1φ2φ3φ400000000φ1φ2φ3φ4u11u21u31u41u12u22u32u42+ξη00ξηλ1λ2.The element pressure is interpolated as follows:(65)p=p1ζ1,where ζ1=1.We discretize equation (48) on the unit Ωj using the interpolation functions of u and p. Since v=v1,v2T∈Xhandq∈Mh is arbitrary, we can choose the following:(1)
v=φl0,l=1,2,3,4, and v=ξη0.(66)ν∇u,∇vj=ν∬Ωj∇u1⋅∇v1+∇u2⋅∇v2dΩ=ν∑j=14uj1∬Ωj∇φj∇φldΩ+νλ1∬Ωj∇ξη∇φldΩv1=φl,l=1,2,3,4,ν∑j=14uj1∬Ωj∇φj∇ξηdΩ+νλ1∬Ωj∇ξη∇ξηdΩv1=ξη,∇⋅v,pj=∬Ωjp∂v1∂x1+∂v2∂x2dΩ=p1∬Ωjζ1∂φl∂x1dΩv1=φl,l=1,2,3,4,p1∬Ωjζ1∂ξη∂x1dΩv1=ξη,f,vj=∬Ωjf1v1+f2v2dΩ=∬Ωjf1φldΩv1=φl,l=1,2,3,4,∬Ωjf1ξηdΩv1=ξη.(2)
v=0φl,l=1,2,3,4, and v=0ξη.(67)ν∇u,∇vj=ν∑j=14uj2∬Ωj∇φj∇φldΩ+νλ1∬Ωj∇ξη∇φldΩv2=φl,l=1,2,3,4,ν∑j=14uj2∬Ωj∇φj∇ξηdΩ+νλ1∬Ωj∇ξη∇ξηdΩv2=ξη,∇⋅v,pj=p1∬Ωjζ1∂φl∂x2dΩv2=φl,l=1,2,3,4,p1∬Ωjζ1∂ξη∂x2dΩv2=ξη,f,vj=∬Ωjf2φldΩv2=φl,l=1,2,3,4,∬Ωjf2ξηdΩv2=ξη.(3)
q=ζ1.(68)∇⋅u,qj=∑j=14uj1∬Ωj∂φj∂x1ζ1dΩ+λ1∬Ωj∂ξη∂x1ζ1dΩ+∑j=14uj2∬Ωj∂φj∂x2ζ1dΩ+λ2∬Ωj∂ξη∂x2ζ1dΩ.By synthesizing (1)–(3) and making some adjustments, the element discrete equations are obtained:(69)∑i=14Kji11ui1+∑i=14Kji12ui2+Kj113λ1+Kj214λ2+Kj115p1=Rj1∑i=14Kji21ui1+∑i=14Kji22ui2+Kj123λ1+Kj224λ2+Kj125p1=Rj2∑i=14K1i31ui1+∑i=14K1i32ui2+K1133λ1+K1234λ2+K1135p1=Rλ11∑i=14K2i41ui1+∑i=14K2i42ui2+K2143λ1+K2244λ2+K2145p1=Rλ22∑i=14K1i51ui1+∑i=14K1i52ui2+K1153λ1+K1254λ2+K1155p1=R13j=1,2,⋯,4.Let(70)δji=1j=i,0j≠i,where(71)Kjiτω=δτων∬Ωj∇φi∇φjdΩτ,ω=1,2Kj113=K1j31=Kj224=K2j42=ν∬Ωj∇ξη∇φjdΩKj214=K2j41=Kj123=K1j32=K1234=K2143=0K1133=K2244=ν∬Ωj∇ξη∇ξηdΩKj1k5=−∬Ωjζ1∂φj∂xkdΩk=1,2K1135=−∬Ωjζ1∂ξη∂x1dΩK2145=−∬Ωjζ1∂ξη∂x2dΩK1i5k=∬Ωj∂φi∂xkζ1dΩk=1,2K1153=∬Ωj∂ξη∂x1ζ1dΩK1254=∬Ωj∂ξη∂x2ζ1dΩK1155=0Rjk=∬ΩjfkφjdΩk=1,2R13=0.Since the velocityuλ in the element does not participate in the distribution of known loads, the equivalent load for λ is zero, i.e., Rλ11=Rλ22=0.The algebraic equations are obtained as follows:(72)K1111⋯K14110⋯0K11130K1115⋮⋱⋮⋮⋱⋮⋮⋮⋮K4111⋯K44110⋯0K41130K41150⋯0K1122⋯K14220K1224K1125⋮⋱⋮⋮⋱⋮⋮⋮⋮0⋯0K4122⋯K44220K4224K4125K1131⋯K14310⋯0K11330K11350⋯0K2142⋯K24420K2244K2145K1151⋯K1451K1152⋯K1452K1153K12540u11⋮u41u12⋮u42λ1λ2p1=R11⋮R41R12⋮R42000.Abbreviations are as follows:(73)K110K130K150K220K24K25K310K330K350K420K44K45K51K52K53K540U1U2λ1λ2P=R1R2000.For the convenience of description, we write the stiffness equation of the mixed element as follows:(74)KaaKaλTKpaTKaλKλλKpλTKpaKpλ0UλP=R00.By eliminating the internal degrees of freedom in the upper first formula and third formula, the condensed element equations are obtained:(75)Kaa−KaλTKλλ−1KaλKpaT−KaλTKλλ−1KpλTKpa−KpλKλλ−1Kaλ−KpλKλλ−1KpλTUP=R0.The coefficient matrix in the equation is the element stiffness matrix with an additional internal velocity term, and the right-hand term is the load column array.


### 3.2. Step 2: Synthesis of Total Stiffness

Let the dimension of the finite element spaceXh and Mh be, respectively, NG1 and NG2, and the basis functions of these two spaces are, respectively, ϕ1,ϕ2,…,ϕNG1 and ψ1,ψ2,…,ψNG2. Then, the solution of the equation can be expressed as(76)uhj=∑i=1NG1uhijϕi,j=1,2ph=∑i=1NG2phiψi.The element stiffness matrix and the load array obtained after condensation are synthesized, respectively. In order to guarantee the uniqueness of the solution ofph, the limit of the pressure Ph is added as follows:(77)∑j=1NG2phj∬ΩψjdΩ=0.The following algebraic equations are synthesized:(78)A11A12A13A21A22A23A31A32A3300A43Uh1Uh2Ph=F1F200.Since the basis functions are piecewise polynomials defined on quasi-uniform regular partitionτh, the coefficient matrix in the system of equation (78) is a large sparse matrix. We use the column principal Gauss elimination method to solve the system. For ease of calculation, we ordered α=Uh1,Uh2,PhT. It comes down to just finding a vector α and setting the coefficient matrix to A1 and B1 for the right-hand vector. Since A1 are 2NG1+NG2+1 rows and 2NG1+NG2 columns, the column principal Gauss elimination method requires that the coefficient matrix must be a square matrix. There are two ways to convert a matrix into a square matrix. The first is to use the least squares to convert A1 into a square matrix of 2NG1+NG2; that is, A1T×A1α=A1T×B1, and the second is to add the last row of A1 to its last column, set the diagonal element to 1, and add a variable, i.e.,(79)A11A12A130A21A22A230A31A32A33A4300A431Uh1Uh2Phx=F1F200.By transformingA1 into 2NG1+NG2+1 matrix, we can obtain the solution of the Stokes problem by using the column principal Gauss elimination method.(1)
Types of finite element:(i)
Douglas rotated element shape function:For the velocityu, we use the four-node Douglas rotated rectangular element to divide the region Ω into several combinations of rectangular elements, each element as shown in Figure 1. Let ai=xi,yi,i=1,2,3,4; then, the centroid coordinates xc,yc are(80)xc=14x1+x2+x3+x4,yc=14y1+y2+y3+y4.The standard element is a square region of spaceξ,η with sides of 2: D:−1≤ξ≤1and−1≤η≤1. The centroid coordinates are at the origin 0,0, and its degree of freedom is the value of four points on a square, which are distributed as shown in Figure 1.Coordinate transformation:(81)ξ=x−xca,η=y−ycb.Convert the rectangular elementΩj in the x,y-plane to the standard element K^ in the ξ,η-plane, where 2a=|x2−x4|and2b=|y3−y1|. Its Jacobi determinant is(82)|J|=|∂x,y∂ξ,η|=ab.Let the shape function of each node ofa^1=0,−1,a^2=1,0,a^3=0,1,anda^4=−1,0 be(83)φ1ξ,η=14−12η+38ξ2−53ξ4−η2−53η4,φ2ξ,η=14+12ξ−38ξ2−53ξ4−η2−53η4,φ3ξ,η=14+12η+38ξ2−53ξ4−η2−53η4,φ4ξ,η=14−12ξ−38ξ2−53ξ4−η2−53η4.The additional bubble function isξη.(ii)
Piecewise constant shape function:For the pressurep, we use a piecewise constant rectangular element to divide the region Ω into several combinations of rectangular elements, each of which is shown in Figure 2. Let a1=xc,yc also be the centroid coordinates. The standard rectangular element is a square region in space ξ,η with a side length of 2, and let a^1=0,0 also be the centroid coordinates at the origin, as shown in Figure 2.The coordinate transformation from rectangular elementΩj to standard element K^ is the same as for velocity u. The node a^1=0,0 corresponds to a shape function(84)ζ1ξ,η=1.(2)
Automatic mesh generation of the finite element method:(i)
Automatic mesh generation of the finite element method for velocityu:In this section, we mainly discuss the automatic generation of 4-node rectangular elements, as well as the node local encoding and the whole encoding comparison table arrayII4,LEE and the node’s actual coordinates’ array XY2,NG1. Let the region Ω be a rectangle, evenly divided by n vertical lines and m horizontal lines; then, the region Ω is divided into n×m quadrilateral elements, as shown in Figure 3.Figure3 represents the overall number and the local number of velocity nodes on a small element. The code of the element and the code of the node are arranged from left to right and from bottom to top, so there are n×m elements and nm+1+n+1m nodes in the whole area Ω. Use i,j to denote the number of vertical and horizontal lines, respectively. Then, i changes from 0 to n,andj changes from 0 to m. The coordinates of i,j for any node are marked xi,yj. Assuming that the middle node of the opposite side of each element is connected by a straight line, there are 2n+1 vertical lines and 2m+1 horizontal lines in the Ω region. So, when i changes from 0 to nandj changes from 0 to 2m, the coordinates xi,yj can be represented by left-right curved margins XLONG and up-down curved margins YLONG. When j is even,(85)xi=XLONGn∗2i+12,yj=YLONGm∗j2.Whenj is odd,(86)xi=XLONGn∗i,yj=YLONGm∗j2.The element whose verticali element band intersects the horizontal j element band is(87)LE=j∗m+i.The number of previous nodes in thej element band is(88)Nj=j∗n+j∗n+1=2n+1∗j.Thus, the overall code for the four nodes of elementLE is(89)II0,LE=Nj+i;II1,LE=Nj+n+i+1II2,LE=Nj+2n+1+i;II3,LE=Nj+n+i.(ii)
Automatic mesh generation of the finite element method for pressurep:In this section, we focus on the automatic generation of 1-node rectangular elements, as well as the node’s actual coordinate arrayXYC2,NG2. If the region Ω is a rectangle, divided evenly by an n vertical line and an m horizontal line, then the region Ω is divided into n×m quadrilateral elements, as shown in Figure 4.Figure4 represents the overall number and the local number of pressure nodes on a small element. The code of the element and the code of the node are arranged from left to right and from bottom to top, so there are n×m elements and nodes in the whole region Ω. Using i,j to denote the number of vertical and horizontal lines, respectively, then i changes from 0 to n,j changes from 0 to m, and each node i,j is marked with the coordinates xi,yj which can be represented by the left-right curve space XLONG and the up-down curve space YLONG.(90)xi=XLONGn∗2i+12,yj=YLONGm∗2j+12.Figure 1
Transformation from the rectangular element to the standard element for the Douglas rotation rectangular element. (a) Rectangular elementΩi. (b) Standard element K^.
(a)(b)Figure 2
Transformation from the rectangular element to the standard element for the piecewise constant rectangular element. (a) Rectangular elementΩi. (b) Standard element K^.
(a)(b)Figure 3
Overall number and local number of velocity nodes. (a) Overall number of nodes. (b) Local number of nodes.
(a)(b)Figure 4
Overall number and local number of pressure nodes. (a) Overall number of nodes. (b) Local number of nodes.
(a)(b)We have given all the constructions of the matrix in the specific form through the overall coding and local coding and the relationship between the base function and the shape function. We can convert all the integrals into standard units, and then we can do all kinds of integrals on standard units.We need to deal with the integral. When we do numerical integral, we usually do not use the formula to calculate the integral, but use the higher algebraic precision Gauss-type numerical integral formula. It is not only convenient to write the standard procedure of the integral but also, according to the construction characteristic of the basis function, it is mainly the integral of the polynomial. If the integral is calculated by the formula, it will cause too much rounding error because of the complicated calculation method. Gauss-type integral is to select some special integral pointξi and weight Hi so that the calculation formula(91)∫−11fξdξ=∑i=1nHifξiis accurate for all 2n−1th polynomials.We have the same formula for double integrals:(92)I=∫−11∫−11fξ,ηdξdη,and then ∫−11fξ,ηdξ=∑j=1nHifξi,η=ϕη.(93)I=∫−11ϕηdη=∑j=1nHjϕηj=∑j=1nHj∑i=1nHifξi,ηj.So,(94)I=∫−11∫−11fξ,ηdξdη=∑i,j=1nHiHjfξi,ηj.With the above analysis, we convert the integrals into double integrals on the standard element and use the 9-node Gauss integral formula; then,(95)∫−11∫−11fξ,ηdξdη=∑i,j=1nHiHjfξi,ηj=∑i,j=1nHijfξi,ηj.The coordinates of Gauss pointξi,ηj and corresponding weights are(96)ξi,ηj=HXi,HXj1≤i,j≤3Hij=GHi∗GHj1≤i,j≤3,where(97)HXi=−0.774597,0,0.7745971≤i,j≤3GHi=0.555556,0.888889,0.5555561≤i,j≤3.


## 3.1. Step 1: Condensation of Internal Degrees of Freedom

Let the four edges of the standard cellK^=−1,12 be, respectively, e^1,e^2,e^3,ande^4; the corresponding midpoints are a^1=0,−1,a^2=1,0,a^3=0,1,anda^4=−1,0, and the shape functions for each node of the Douglas rotation element obtained are, respectively,(60)φ1ξ,η=14−12η+38ξ2−53ξ4−η2−53η4φ2ξ,η=14+12ξ−38ξ2−53ξ4−η2−53η4φ3ξ,η=14+12η+38ξ2−53ξ4−η2−53η4φ4ξ,η=14−12ξ−38ξ2−53ξ4−η2−53η4.It is easy to verify(61)∫e^iφjds^=δij|e^i|,i,j=1,2,3,4,where δij is the Kronecker symbol and e^i is the length of the edge of e^i.The element speed interpolation with additional node-free speed terms can be expressed as follows:(62)u1=∑i=14φiui1+λ1ξη,u2=∑i=14φiui2+λ2ξη,where λ1andλ2 are the undetermined coefficient in the element’s internal definition, known as the internal degree of freedom.The expression is expressed in the matrix form:(63)u=ua+uλ,i.e.,(64)u1u2=φ1φ2φ3φ400000000φ1φ2φ3φ4u11u21u31u41u12u22u32u42+ξη00ξηλ1λ2.The element pressure is interpolated as follows:(65)p=p1ζ1,where ζ1=1.We discretize equation (48) on the unit Ωj using the interpolation functions of u and p. Since v=v1,v2T∈Xhandq∈Mh is arbitrary, we can choose the following:(1)
v=φl0,l=1,2,3,4, and v=ξη0.(66)ν∇u,∇vj=ν∬Ωj∇u1⋅∇v1+∇u2⋅∇v2dΩ=ν∑j=14uj1∬Ωj∇φj∇φldΩ+νλ1∬Ωj∇ξη∇φldΩv1=φl,l=1,2,3,4,ν∑j=14uj1∬Ωj∇φj∇ξηdΩ+νλ1∬Ωj∇ξη∇ξηdΩv1=ξη,∇⋅v,pj=∬Ωjp∂v1∂x1+∂v2∂x2dΩ=p1∬Ωjζ1∂φl∂x1dΩv1=φl,l=1,2,3,4,p1∬Ωjζ1∂ξη∂x1dΩv1=ξη,f,vj=∬Ωjf1v1+f2v2dΩ=∬Ωjf1φldΩv1=φl,l=1,2,3,4,∬Ωjf1ξηdΩv1=ξη.(2)
v=0φl,l=1,2,3,4, and v=0ξη.(67)ν∇u,∇vj=ν∑j=14uj2∬Ωj∇φj∇φldΩ+νλ1∬Ωj∇ξη∇φldΩv2=φl,l=1,2,3,4,ν∑j=14uj2∬Ωj∇φj∇ξηdΩ+νλ1∬Ωj∇ξη∇ξηdΩv2=ξη,∇⋅v,pj=p1∬Ωjζ1∂φl∂x2dΩv2=φl,l=1,2,3,4,p1∬Ωjζ1∂ξη∂x2dΩv2=ξη,f,vj=∬Ωjf2φldΩv2=φl,l=1,2,3,4,∬Ωjf2ξηdΩv2=ξη.(3)
q=ζ1.(68)∇⋅u,qj=∑j=14uj1∬Ωj∂φj∂x1ζ1dΩ+λ1∬Ωj∂ξη∂x1ζ1dΩ+∑j=14uj2∬Ωj∂φj∂x2ζ1dΩ+λ2∬Ωj∂ξη∂x2ζ1dΩ.By synthesizing (1)–(3) and making some adjustments, the element discrete equations are obtained:(69)∑i=14Kji11ui1+∑i=14Kji12ui2+Kj113λ1+Kj214λ2+Kj115p1=Rj1∑i=14Kji21ui1+∑i=14Kji22ui2+Kj123λ1+Kj224λ2+Kj125p1=Rj2∑i=14K1i31ui1+∑i=14K1i32ui2+K1133λ1+K1234λ2+K1135p1=Rλ11∑i=14K2i41ui1+∑i=14K2i42ui2+K2143λ1+K2244λ2+K2145p1=Rλ22∑i=14K1i51ui1+∑i=14K1i52ui2+K1153λ1+K1254λ2+K1155p1=R13j=1,2,⋯,4.Let(70)δji=1j=i,0j≠i,where(71)Kjiτω=δτων∬Ωj∇φi∇φjdΩτ,ω=1,2Kj113=K1j31=Kj224=K2j42=ν∬Ωj∇ξη∇φjdΩKj214=K2j41=Kj123=K1j32=K1234=K2143=0K1133=K2244=ν∬Ωj∇ξη∇ξηdΩKj1k5=−∬Ωjζ1∂φj∂xkdΩk=1,2K1135=−∬Ωjζ1∂ξη∂x1dΩK2145=−∬Ωjζ1∂ξη∂x2dΩK1i5k=∬Ωj∂φi∂xkζ1dΩk=1,2K1153=∬Ωj∂ξη∂x1ζ1dΩK1254=∬Ωj∂ξη∂x2ζ1dΩK1155=0Rjk=∬ΩjfkφjdΩk=1,2R13=0.Since the velocityuλ in the element does not participate in the distribution of known loads, the equivalent load for λ is zero, i.e., Rλ11=Rλ22=0.The algebraic equations are obtained as follows:(72)K1111⋯K14110⋯0K11130K1115⋮⋱⋮⋮⋱⋮⋮⋮⋮K4111⋯K44110⋯0K41130K41150⋯0K1122⋯K14220K1224K1125⋮⋱⋮⋮⋱⋮⋮⋮⋮0⋯0K4122⋯K44220K4224K4125K1131⋯K14310⋯0K11330K11350⋯0K2142⋯K24420K2244K2145K1151⋯K1451K1152⋯K1452K1153K12540u11⋮u41u12⋮u42λ1λ2p1=R11⋮R41R12⋮R42000.Abbreviations are as follows:(73)K110K130K150K220K24K25K310K330K350K420K44K45K51K52K53K540U1U2λ1λ2P=R1R2000.For the convenience of description, we write the stiffness equation of the mixed element as follows:(74)KaaKaλTKpaTKaλKλλKpλTKpaKpλ0UλP=R00.By eliminating the internal degrees of freedom in the upper first formula and third formula, the condensed element equations are obtained:(75)Kaa−KaλTKλλ−1KaλKpaT−KaλTKλλ−1KpλTKpa−KpλKλλ−1Kaλ−KpλKλλ−1KpλTUP=R0.The coefficient matrix in the equation is the element stiffness matrix with an additional internal velocity term, and the right-hand term is the load column array.


## 3.2. Step 2: Synthesis of Total Stiffness

Let the dimension of the finite element spaceXh and Mh be, respectively, NG1 and NG2, and the basis functions of these two spaces are, respectively, ϕ1,ϕ2,…,ϕNG1 and ψ1,ψ2,…,ψNG2. Then, the solution of the equation can be expressed as(76)uhj=∑i=1NG1uhijϕi,j=1,2ph=∑i=1NG2phiψi.The element stiffness matrix and the load array obtained after condensation are synthesized, respectively. In order to guarantee the uniqueness of the solution ofph, the limit of the pressure Ph is added as follows:(77)∑j=1NG2phj∬ΩψjdΩ=0.The following algebraic equations are synthesized:(78)A11A12A13A21A22A23A31A32A3300A43Uh1Uh2Ph=F1F200.Since the basis functions are piecewise polynomials defined on quasi-uniform regular partitionτh, the coefficient matrix in the system of equation (78) is a large sparse matrix. We use the column principal Gauss elimination method to solve the system. For ease of calculation, we ordered α=Uh1,Uh2,PhT. It comes down to just finding a vector α and setting the coefficient matrix to A1 and B1 for the right-hand vector. Since A1 are 2NG1+NG2+1 rows and 2NG1+NG2 columns, the column principal Gauss elimination method requires that the coefficient matrix must be a square matrix. There are two ways to convert a matrix into a square matrix. The first is to use the least squares to convert A1 into a square matrix of 2NG1+NG2; that is, A1T×A1α=A1T×B1, and the second is to add the last row of A1 to its last column, set the diagonal element to 1, and add a variable, i.e.,(79)A11A12A130A21A22A230A31A32A33A4300A431Uh1Uh2Phx=F1F200.By transformingA1 into 2NG1+NG2+1 matrix, we can obtain the solution of the Stokes problem by using the column principal Gauss elimination method.(1)
Types of finite element:(i)
Douglas rotated element shape function:For the velocityu, we use the four-node Douglas rotated rectangular element to divide the region Ω into several combinations of rectangular elements, each element as shown in Figure 1. Let ai=xi,yi,i=1,2,3,4; then, the centroid coordinates xc,yc are(80)xc=14x1+x2+x3+x4,yc=14y1+y2+y3+y4.The standard element is a square region of spaceξ,η with sides of 2: D:−1≤ξ≤1and−1≤η≤1. The centroid coordinates are at the origin 0,0, and its degree of freedom is the value of four points on a square, which are distributed as shown in Figure 1.Coordinate transformation:(81)ξ=x−xca,η=y−ycb.Convert the rectangular elementΩj in the x,y-plane to the standard element K^ in the ξ,η-plane, where 2a=|x2−x4|and2b=|y3−y1|. Its Jacobi determinant is(82)|J|=|∂x,y∂ξ,η|=ab.Let the shape function of each node ofa^1=0,−1,a^2=1,0,a^3=0,1,anda^4=−1,0 be(83)φ1ξ,η=14−12η+38ξ2−53ξ4−η2−53η4,φ2ξ,η=14+12ξ−38ξ2−53ξ4−η2−53η4,φ3ξ,η=14+12η+38ξ2−53ξ4−η2−53η4,φ4ξ,η=14−12ξ−38ξ2−53ξ4−η2−53η4.The additional bubble function isξη.(ii)
Piecewise constant shape function:For the pressurep, we use a piecewise constant rectangular element to divide the region Ω into several combinations of rectangular elements, each of which is shown in Figure 2. Let a1=xc,yc also be the centroid coordinates. The standard rectangular element is a square region in space ξ,η with a side length of 2, and let a^1=0,0 also be the centroid coordinates at the origin, as shown in Figure 2.The coordinate transformation from rectangular elementΩj to standard element K^ is the same as for velocity u. The node a^1=0,0 corresponds to a shape function(84)ζ1ξ,η=1.(2)
Automatic mesh generation of the finite element method:(i)
Automatic mesh generation of the finite element method for velocityu:In this section, we mainly discuss the automatic generation of 4-node rectangular elements, as well as the node local encoding and the whole encoding comparison table arrayII4,LEE and the node’s actual coordinates’ array XY2,NG1. Let the region Ω be a rectangle, evenly divided by n vertical lines and m horizontal lines; then, the region Ω is divided into n×m quadrilateral elements, as shown in Figure 3.Figure3 represents the overall number and the local number of velocity nodes on a small element. The code of the element and the code of the node are arranged from left to right and from bottom to top, so there are n×m elements and nm+1+n+1m nodes in the whole area Ω. Use i,j to denote the number of vertical and horizontal lines, respectively. Then, i changes from 0 to n,andj changes from 0 to m. The coordinates of i,j for any node are marked xi,yj. Assuming that the middle node of the opposite side of each element is connected by a straight line, there are 2n+1 vertical lines and 2m+1 horizontal lines in the Ω region. So, when i changes from 0 to nandj changes from 0 to 2m, the coordinates xi,yj can be represented by left-right curved margins XLONG and up-down curved margins YLONG. When j is even,(85)xi=XLONGn∗2i+12,yj=YLONGm∗j2.Whenj is odd,(86)xi=XLONGn∗i,yj=YLONGm∗j2.The element whose verticali element band intersects the horizontal j element band is(87)LE=j∗m+i.The number of previous nodes in thej element band is(88)Nj=j∗n+j∗n+1=2n+1∗j.Thus, the overall code for the four nodes of elementLE is(89)II0,LE=Nj+i;II1,LE=Nj+n+i+1II2,LE=Nj+2n+1+i;II3,LE=Nj+n+i.(ii)
Automatic mesh generation of the finite element method for pressurep:In this section, we focus on the automatic generation of 1-node rectangular elements, as well as the node’s actual coordinate arrayXYC2,NG2. If the region Ω is a rectangle, divided evenly by an n vertical line and an m horizontal line, then the region Ω is divided into n×m quadrilateral elements, as shown in Figure 4.Figure4 represents the overall number and the local number of pressure nodes on a small element. The code of the element and the code of the node are arranged from left to right and from bottom to top, so there are n×m elements and nodes in the whole region Ω. Using i,j to denote the number of vertical and horizontal lines, respectively, then i changes from 0 to n,j changes from 0 to m, and each node i,j is marked with the coordinates xi,yj which can be represented by the left-right curve space XLONG and the up-down curve space YLONG.(90)xi=XLONGn∗2i+12,yj=YLONGm∗2j+12.Figure 1
Transformation from the rectangular element to the standard element for the Douglas rotation rectangular element. (a) Rectangular elementΩi. (b) Standard element K^.
(a)(b)Figure 2
Transformation from the rectangular element to the standard element for the piecewise constant rectangular element. (a) Rectangular elementΩi. (b) Standard element K^.
(a)(b)Figure 3
Overall number and local number of velocity nodes. (a) Overall number of nodes. (b) Local number of nodes.
(a)(b)Figure 4
Overall number and local number of pressure nodes. (a) Overall number of nodes. (b) Local number of nodes.
(a)(b)We have given all the constructions of the matrix in the specific form through the overall coding and local coding and the relationship between the base function and the shape function. We can convert all the integrals into standard units, and then we can do all kinds of integrals on standard units.We need to deal with the integral. When we do numerical integral, we usually do not use the formula to calculate the integral, but use the higher algebraic precision Gauss-type numerical integral formula. It is not only convenient to write the standard procedure of the integral but also, according to the construction characteristic of the basis function, it is mainly the integral of the polynomial. If the integral is calculated by the formula, it will cause too much rounding error because of the complicated calculation method. Gauss-type integral is to select some special integral pointξi and weight Hi so that the calculation formula(91)∫−11fξdξ=∑i=1nHifξiis accurate for all 2n−1th polynomials.We have the same formula for double integrals:(92)I=∫−11∫−11fξ,ηdξdη,and then ∫−11fξ,ηdξ=∑j=1nHifξi,η=ϕη.(93)I=∫−11ϕηdη=∑j=1nHjϕηj=∑j=1nHj∑i=1nHifξi,ηj.So,(94)I=∫−11∫−11fξ,ηdξdη=∑i,j=1nHiHjfξi,ηj.With the above analysis, we convert the integrals into double integrals on the standard element and use the 9-node Gauss integral formula; then,(95)∫−11∫−11fξ,ηdξdη=∑i,j=1nHiHjfξi,ηj=∑i,j=1nHijfξi,ηj.The coordinates of Gauss pointξi,ηj and corresponding weights are(96)ξi,ηj=HXi,HXj1≤i,j≤3Hij=GHi∗GHj1≤i,j≤3,where(97)HXi=−0.774597,0,0.7745971≤i,j≤3GHi=0.555556,0.888889,0.5555561≤i,j≤3.


## 4. A Nonconforming Rotated Finite Element Method for the Navier–Stokes Equations

In the first section of this section, the existence, uniqueness, and stability of the solution of the NRFE approximation problem are discussed. The convergence rate of the NRFE method is given in the second section; there are two steps in solving the Navier–Stokes equations with the NRFE method: dealing with the linear part and dealing with the nonlinear part. The first step has been solved in the last section. Here, we mainly discuss how to solve the nonlinear term by the Newton iteration method.Let⋅,⋅j=⋅,⋅Ωj,1≤j≤J, define the discrete bilinear form(98)ahu,v=ν∑j=1J∇u,∇vj,∀u,v∈X∪Xhdhv,q=∑j=1J∇⋅v,qj,∀v∈Xh,q∈Mh.and norm(99)‖v‖1,h=∑j=1J‖∇v‖0,j21/2=1νahv,v,∀v∈Xh.For∀u,v,w∈X∪Xh, define the discrete trilinear form(100)b1,hu,v,w=∑j=1Ju⋅∇v,wjbhu,v,w=12b1,hu,v,w−b1,hu,w,v,which satisfies the following properties:(i)
(101)b1,hu,v,w+b1,hu,w,v=∑j=1Ju⋅njv,wj−∑j=1J∇⋅uv,wj.(ii)
(102)b1,hu,v,w=bhu,v,w+12∑j=1Ju⋅njv,wj−12∑j=1J∇⋅uv,wj.(iii)
(103)b1,hu,v,w≤cb‖u‖1,h‖v‖1,h‖w‖1,h,bhu,v,w≤cb‖u‖1,h‖v‖1,h‖w‖1,h.(iv)
(104)bhu,v,w=−bhu,w,v.Here and laternj,1≤j≤J, represents the unit outer normal vector of ∂Ωj. See [7, 17] for its proofs.In addition, we introduce a closed subset ofXh:(105)Vh=v∈Xh;dhv,q=0,∀q∈Mh.Then, the nonconforming rotated finite element approximation problem of problem (6) is(106)finduh,ph∈Xh×Mh,soahuh,v−dhv,ph+bhuh,uh,v=f,v,∀v∈Xh,dhuh,q=0,∀q∈Mh.Theorem 7 (see [7]).
Under the condition that Proposition1 is established, if(107)cbν2‖f‖∗<1,‖f‖∗=supv∈Xhf,v‖v‖1,h,then there is, to problem (50), a unique solution uh,ph∈Xh×Mh, and uh∈Vh satisfies(108)‖uh‖1,h≤1ν‖f‖∗.Proposition 4 (see [7]).
Letu,p and uh,ph be the solutions of problems (11) and (106), respectively; then, under the condition of formula (51), there exists a positive number c such that(109)‖u−uh‖1,h≤cinfv∈Xh‖u−v‖1,h+supv∈Vhahu,v+bhu,u,v−f,v‖v‖1,h‖p−ph‖0≤cinfq∈Mh‖p−q‖0+supv∈Xhahu,v−dhv,p+bhu,u,v−f,v‖v‖1,h+infv∈Xh‖u−v‖1,h+supv∈Vhahu,v+bhu,u,v−f,v‖v‖1,h.Proposition 5 (see [7]).
Letf,gj=∫∂Ωjfgds, for ∀ϕ,w∈X∪Xh, satisfy the following properties:(110)|∑j=1J∂w∂nj,ϕj|≤ch‖w‖2‖ϕ‖1,h,∀w∈DA|∑j=1Jw⋅njv,ϕj|≤ch‖w‖1,h‖v‖2‖ϕ‖1,h,∀v∈H2Ω2|∑j=1Jq,ϕ⋅njj|≤ch‖q‖1‖ϕ‖1,h,∀q∈H1Ω.Theorem 8 (see [7]).
Letu,p∈H2Ω2×H1Ω and uh,ph∈Xh×Mh be the solutions of problems (11) and (106), respectively; then,(111)‖u−uh‖1,h+‖p−ph‖0≤ch‖u‖2+‖p‖1.Using the dual argument method, we estimate the error ofL2.Consider the linear duality problem:(112)−νΔψ+∇χ+u⋅∇ψ+∇u⋅ψ=u−uh∀x∈Ω,∇⋅ψ=0∀x∈Ω,ψ=0∀x∈Γ.Its variational form is(113)findψ,χ∈X×M,soav,ψ−dv,χ+u⋅∇v,ψ+v⋅∇u,ψ=u−uh,v,∀v∈X,dψ,q=0,∀q∈M.Ifu is the nonsingular solution of problem (1), then problem (113) has a unique solution (see [16]).In order to establish an estimate of theL2 error of velocity, the Aubin–Nitsche argumentation method is used here. Let problem (112) be H2 regular, i.e.,(114)‖ψ‖2+‖χ‖1≤c‖u−uh‖0.Letψ,χ be the solution of problem (113), and ψh,χh∈Xh×Mh satisfies(115)‖ψ−ψh‖1,h+‖χ−χh‖0≤ch‖ψ‖2+‖χ‖1.Theorem 9 (see [7]).
Letu,p∈H2Ω2×H1Ω and uh,ph∈Xh×Mh be the solutions of problems (11) and (50), respectively; if problem (112) is H2 regular, then(116)‖u−uh‖0≤ch2‖u‖2+‖p‖1.There are two steps in solving the Navier–Stokes equations using the NRFE method: dealing with the linear part and dealing with the nonlinear part. The method of dealing with the linear term is the same as that in Section3. When dealing with the nonlinear term, the most common method is the Newton iterative method, and the Newton iterative method requires higher initial value, so we choose the solution of the Stokes problem as the initial value. Let the initial value of the iteration be uh0,ph0 and uhn,phn be the nth computing solution; then, the iteration format is as follows:(117)ahuhn+1,v−dhv,phn+1+bhuhn+1,uhn,v+bhuhn,uhn+1,v=f,v+bhuhn,uhn,v,∀v∈Xh,dhuhn+1,q=0,∀q∈Mh.In this part, we only need to find the matrix obtained from the nonlinear part and then add it to the matrix obtained from the linear part (i.e., the coefficient matrix of the Stokes equationsA1 and the right-hand term B1). Then, the solution of the Navier–Stokes equations can be obtained by the iterative method. In the following, we come to explain the concrete process of discretization.(1)
Usingv=ϕl0,l=1,…,NG1 to discrete equation (117), the iterative partial result is(118)bhuhn+1,uhn,v=12∑i=1NG1∑j=1NG1un+1i1unj1∬Ωϕi∂ϕj∂x1ϕl−ϕi∂ϕl∂x1ϕjdΩ+un+1i2unj1∬Ωϕi∂ϕj∂x2ϕl−ϕi∂ϕl∂x2ϕjdΩ,bhuhn,uhn+1,v=12∑i=1NG1∑j=1NG1unj1un+1i1∬Ωϕj∂ϕi∂x1ϕl−ϕj∂ϕl∂x1ϕidΩ+unj2un+1i1∬Ωϕj∂ϕi∂x2ϕl−ϕj∂ϕl∂x2ϕidΩ,bhuhn,uhn,v=12∑i=1NG1∑j=1NG1uni1unj1∬Ωϕi∂ϕj∂x1ϕl−ϕi∂ϕl∂x1ϕjdΩ+uni2unj1∬Ωϕi∂ϕj∂x2ϕl−ϕi∂ϕl∂x2ϕjdΩ.(2)
Usingv=0ϕl,l=1,…,NG1 to discrete equation (117), the iterative partial result is(119)bhuhn+1,uhn,v=12∑i=1NG1∑j=1NG1un+1i1unj2∬Ωϕi∂ϕj∂x1ϕl−ϕi∂ϕl∂x1ϕjdΩ+un+1i2unj2∬Ωϕi∂ϕj∂x2ϕl−ϕi∂ϕl∂x2ϕjdΩ,bhuhn,uhn+1,v=12∑i=1NG1∑j=1NG1unj1un+1i2∬Ωϕj∂ϕi∂x1ϕl−ϕj∂ϕl∂x1ϕidΩ+unj2un+1i2∬Ωϕj∂ϕi∂x2ϕl−ϕj∂ϕl∂x2ϕidΩ,bhuhn,uhn,v=12∑i=1NG1∑j=1NG1uni1unj2∬Ωϕi∂ϕj∂x1ϕl−ϕi∂ϕl∂x1ϕjdΩ+uni2unj2∬Ωϕi∂ϕj∂x2ϕl−ϕi∂ϕl∂x2ϕjdΩ.After discretization, the resulting iteration is as follows:(120)A2=A11A120A21A220000000B2=F1F200,where(121)Aτω=aliNG1×NG1,Fτ=blNG1ali=12∑j=1NG1unjσ∬Ωϕi∂ϕj∂xσϕl−ϕi∂ϕl∂xσϕjdΩ+12∑j=1NG1unjσ∬Ωϕj∂ϕi∂xσϕl−ϕj∂ϕl∂xσϕidΩτ=ω12∑j=1NG1unjγ∬Ωϕj∂ϕi∂xγϕl−ϕj∂ϕl∂xγϕidΩ+12∑j=1NG1unjτ∬Ωϕi∂ϕj∂xγϕl−ϕi∂ϕl∂xγϕjdΩτ≠ωbl=12∑j=1NG1∑i=1NG1unjσuniσ∬Ωϕj∂ϕi∂xσϕl−ϕj∂ϕl∂xσϕidΩ+unjγuniσ∬Ωϕj∂ϕi∂xγϕl−ϕj∂ϕl∂xγϕidΩ.In the formula above,τ,ω,σ,γ=1,2, and σ≠γ.In this way, we can construct a finite element system of equations,A1+A2αm+1=B1+B2, and solve it iteratively. Until ‖αm+1−αm‖ is less than given ε, we can get the solution uh,ph of equation (64).


## 5. A Two-Level Nonconforming Rotated Finite Element Method for the Navier–Stokes Equations

Because of the nonlinear and incompressible conditions of the Navier–Stokes equations, it is very difficult to solve them numerically, and they consume a lot of CPU time. To overcome these difficulties, many methods have been proposed, including the two-level method. The aim of the two-level method is to obtain the approximate solution of the nonlinear equation in less time and to maintain the optimal convergence speed. The concrete calculation is to solve a nonlinear problem on the coarse grid and a linear problem on the fine grid. The theoretical and numerical experiments in this paper show that, compared with the NRFE method, the TNRFE method can save a lot of time when the convergence rate reaches the same order, which shows that it is an effective algorithm. In this section, we prove the stability and convergence rate of the TNRFE method.We choose a coarse gridH and a fine grid h with H≫0,H,h to represent the coarse and fine grid scales of the velocity and pressure spaces and then establish a nonconforming finite element space XH×MH, Xh×Mh, and XH×MH⊂Xh×Mh. Similarly, they satisfy all the properties mentioned in the previous section.The Newton two-level nonconforming finite element method is divided into two steps:Step 1: solving a nonlinear Navier–Stokes problem on a coarse grid:(122)finduH,pH∈XH×MH,soahuH,v−dhv,pH+dhuH,q+bhuH,uH,v=f,v,∀v,q∈XH×MH.Step 2: solving a linearized Navier–Stokes problem on a fine grid:(123)finduh,ph∈Xh×Mh,soahuh,v−dhv,ph+dhuh,q+bhuh,uH,v+bhuH,uh,v=f,v+bhuH,uH,v,∀v,q∈Xh×Mh.Define a generalized bilinear form onX×M×X×M:(124)ℬu,p;v,q=au,v−dv,p+du,q.So, the variational form of problem (1) is as follows: find u,p∈X×M such that(125)ℬu,p;v,q+bu,u,v=f,v,∀v,q∈X×M.For givenu,p∈X×M, we define its Stokes projection w,π∈Xh×Mh as(126)ℬhw,π;v,q=ℬu,p;v,q,∀v,q∈Xh×Mh.Then, nonconforming rotated finite element approximation problem (124) is as follows: find uh,ph∈Xh×Mh so that(127)ℬhuh,ph;v,q+bhuh,uh,v=f,v,∀v,q∈Xh×Mh.Proposition 6.
(i)
There is a constantc>0 such that(128)|ℬu,p;v,q|≤c‖∇u‖0+‖p‖0‖∇v‖0+‖q‖0,∀u,p,v,q∈X×M.(ii)
There exists a constantα>0 that depends on h:(129)α‖uh‖1,h+‖ph‖0≤supvh,qh∈Xh×Mhℬhuh,ph;vh,qh‖vh‖1,h+‖qh‖0,∀uh,ph∈Xh×Mh.Proof .
((70) proof (see [4, 16–18])). The following is to prove (129). Let uh,ph be a nonzero element in Xh×Mh. Let qh∈Mh be determined by the problem dhψh,ph=dhψh,qh,∀ψh∈Xh, so(130)ℬhuh,ph;uh,qh=ν‖uh‖1,h2−dhuh,ph+dhuh,qh=ν‖uh‖1,h2.
Since there exists the constantδ>0 for all ph∈Mh⊂M, there exists v∈X satisfying ∇⋅v=−ph and δ‖∇v‖0≤‖ph‖0. Hence, derived from upper equations (39) and (40) and the Young inequality,(131)ℬhuh,ph;Πhv,0≥−ν‖uh‖1,h‖Πhv‖1,h−dhΠhv,ph≥−cν‖uh‖1,h‖∇v‖0−dhΠhv−v,ph+‖ph‖02≥−cν‖uh‖1,h‖∇v‖0+‖ph‖02≥−c2ν22δ‖uh‖1,h2−δ2‖∇v‖02+‖ph‖02≥−c2ν22δ‖uh‖1,h2+12‖ph‖02.
Letwh=uh+δ/c2νΠhv and rh=qh; from the above two formulas, we can get(132)ℬhuh,ph;wh,rh≥ν2‖uh‖1,h2+δ2c2ν‖ph‖02.
Obtained from (50) and (40), the triangle inequality, and the Hölder inequality,(133)‖wh‖1,h+‖rh‖0=‖uh+δc2νΠhv‖1,h+‖qh‖0≤‖uh‖1,h+δc2ν‖Πhv‖1,h+1β0supψh∈Xhdhψh,qh‖ψh‖1,h≤‖uh‖1,h+δcν‖∇v‖0+1β0supψh∈Xhdhψh,ph‖ψh‖1,h≤‖uh‖1,h+1cν‖ph‖0+2β0‖ph‖0.
From the above two formulas, it is easy to get (129).Theorem 10.
Under the conditions of Proposition6 and Theorem 7, the solution uh,ph∈Xh×Mh satisfies(134)‖uh‖1,h+‖ph‖0≤c‖f‖∗.
Here,c is a constant that depends on ν,f,andΩ.Proof.
According to Theorem7 and (122), there exists a unique solution uH,pH∈XH×MH, and it satisfies(135)‖uH‖1,h≤1ν‖f‖∗.
In equation (123), we make v=uhandq=ph; then, from (104), we can derive(136)ν‖uh‖1,h2+bhuh,uH,uh=f,uh+bhuH,uH,uh.
Using (103), we can get(137)ν‖uh‖1,h2≤cb‖uh‖1,h2‖uH‖1,h+‖f‖∗‖uh‖1,h+cb‖uH‖1,h2‖uh‖1,h.
From (135), we can derive(138)ν‖uh‖1,h2≤cbν‖f‖∗‖uh‖1,h2+‖f‖∗‖uh‖1,h+cbν2‖f‖∗2‖uh‖1,h.
After sorting through the equations, we can get(139)ν1−cbν2‖f‖∗‖uh‖1,h≤1+cbν2‖f‖∗‖f‖∗.
From Theorem7, we can get(140)‖uh‖1,h≤c‖f‖∗.
From Proposition6 and (123), we can get(141)α‖uh‖1,h+‖ph‖0≤supv,q∈Xh×Mhℬhuh,ph;v,q‖v‖1,h+‖q‖0≤cb‖uH‖1,h2+‖f‖∗+2cb‖uh‖1,h‖uH‖1,h.
From (135), (140), and (141), we can derive(142)‖uh‖1,h+‖ph‖0≤c‖f‖∗.Theorem 11.
Under the conditions of Proposition6 and Theorem 7, there exists a constant c which has nothing to do with h such that the solution uh,ph of problems (122) and (123) satisfies(143)‖u−uh‖1,h+‖p−ph‖0≤ch+H2.Proof .
Multiply the two ends of the first formula of (1) by v∈Xh and (102), and we can get(144)ahu,v−dhv,p+bhu,u,v−ν∑j=1J∂u∂nj,vj+∑j=1Jp,v⋅njj+12∑j=1Ju⋅nju,vj=f,v.
Supposew,π is the projection of u,p on Xh×Mh, and here, w∈Vh; then, by (125) and (144), we can get(145)ahuh−w,v−dhv,ph−π+dhuh−w,q=ℬhuh−w,ph−π;v,q=bhu,u,v−ν∑j=1J∂u∂nj,vj+∑j=1Jp,v⋅njj+12∑j=1Ju⋅nju,vj+bhuH,uH,v−bhuh,uH,v−bhuH,uh,v.
In (145), we take v=uh−wandq=ph−π. By using (3), (103), and Proposition 5, we can get(146)ν‖uh−w‖1,h2=bhu,u,uh−w−ν∑j=1J∂u∂nj,uh−wj+∑j=1J<p,uh−w⋅nj>j+12∑j=1Ju⋅nju,uh−wj+bhuH,uH,uh−w−bhuh,uH,uh−w−bhuH,uh,uh−w=bhu−w,u,uh−w+bhw,u−w,uh−w−bhuh−w,uH,uh−w−bhuH−w,uh−w,uH−w−ν∑j=1J∂u∂nj,uh−wj+∑j=1Jp,uh−w⋅njj+12∑j=1Ju⋅nju,uh−wj≤cb‖uh−w‖1,h2‖uH‖1,h+cb‖uh−w‖1,h‖u−w‖1,h‖∇u‖0+‖w‖1,h‖u−w‖1,h+‖uH−w‖1,h2+ch‖uh−w‖1,h‖u‖2+‖p‖1+‖∇u‖0‖u‖2.
From Theorems2, 3, 7, and 8 and (135), we can get(147)‖uh−w‖1,h≤ch+H2.
From Proposition6 and (145), we can get(148)‖uh−w‖1,h+‖ph−π‖0≤α−1supv,q∈Xh×Mhℬhuh−w,ph−π;v,q‖v‖1,h+‖q‖0≤c‖∇u‖0+‖uh−w‖1,h+‖w‖1,h‖u−w‖1,h+‖uh−w‖1,h+c‖uh−w‖1,h+‖w−uH‖1,h‖uh−w‖1,h+c‖w−uH‖1,h2+ch‖u‖2+‖p‖1+‖∇u‖0‖u‖2≤ch+H2.
So,(149)‖u−uh‖1,h+‖p−ph‖0≤‖u−w‖1,h+‖w−uh‖1,h+‖p−π‖0+‖π−ph‖0≤ch+H2.
From the above results, when the appropriate scaling relation is chosen, i.e.,h=OH2, we can see that the NRFE method and the Newton TNRFE method have the same convergence rate.


## 6. Computational Example

Here, examples are given to verify the correctness of the theory and to show that the two-level nonconforming finite element method for solving the Navier–Stokes equations is superior to the nonconforming finite element method.Example 1.
Let us take a true solution which satisfies the conditions and set the two-dimensional rectangular regionΩ=0,1×0,1; the true solutions are as follows:(150)ux,y=u1x,y,u2x,y,u1x,y=x2x−12yy−12y−1,u2x,y=−xx−12x−1y2y−12,px,y=x−0.5.
It is proved that the velocity term of this vector function belongs to the nondivergence space, i.e.,∇⋅u=0, and the pressure term satisfies ∫Ωpxdx=0. By substituting ux,y,px,y into Stokes (12), the right-hand term fx,y=f1x,y,f2x,y is obtained. Since the viscosity coefficient has little effect on the accuracy of the equation, we only consider the case when the viscosity coefficient is ν=1. The numerical results of the nonconforming finite element method with different values of h are given in Table 1.
From Table1, for solving the Stokes equations, we can see that the NRFE method converges with h decreasing, and it keeps the convergence rate. When the grid scale is h=1/25, the true solution and the calculation solution of the velocity component u1 are shown in Figure 5.
The contour of the true solution and the calculation solution of the velocity componentu1 is shown in Figure 6.
The true solution and the calculation solution of the velocity componentu2 are shown in Figure 7.
The contour of the true solution and the calculation solution of the velocity componentu2 is shown in Figure 8.
The vector of the true solution and the calculation solution of the velocity is shown in Figure9.
The images of the true solution and the calculation solution of the pressure are shown in Figure10.
The contours of the true solution and the calculated solution of the pressure are shown in Figure11.
Here, examples are given to verify the correctness of the theory and to show that the TNRFE method for solving the Navier–Stokes equations is superior to the NRFE method.Table 1
Convergence study of the NRFE method for the Stokes equations in Example1.
hCPU (s)‖u−uh‖1,h/‖u‖1,h‖u−uh‖0/‖u‖0‖p−ph‖0/‖p‖01/910.052767352940.006833781230.007216578841/16140.049305443400.005959344260.006729089951/252490.046267803790.005242388970.00630336367Figure 5
True solution and calculation solution of velocity componentu1 by the NRFE method for the Stokes equations. (a) True solution of velocity component u1. (b) Calculation solution of velocity component u1.
(a)(b)Figure 6
True solution contour and calculation solution contour of velocity componentu1 by the NRFE method for the Stokes equations. (a) True solution contour of velocity component u1. (b) Calculation solution contour of velocity component u1.
(a)(b)Figure 7
True solution and calculation solution of velocity componentu2 by the NRFE method for the Stokes equations. (a) True solution of velocity component u2. (b) Calculation solution of velocity component u2.
(a)(b)Figure 8
True solution contour and calculation solution contour of velocity componentu2 by the NRFE method for the Stokes equations. (a) True solution contour of velocity component u2. (b) Calculation solution contour of velocity component u2.
(a)(b)Figure 9
True solution vector and calculation solution vector of velocity by the NRFE method for the Stokes equations. (a) True solution vector of velocity. (b) Calculation solution vector of velocity.
(a)(b)Figure 10
True solution and calculation solution of pressure by the NRFE method for the Stokes equations. (a) True solution of pressure. (b) Calculation solution of pressure.
(a)(b)Figure 11
True solution contour and calculation solution contour of pressure by the NRFE method for the Stokes equations. (a) True solution contour of pressure. (b) Calculation solution contour of pressure.
(a)(b)Example 2.
Let us take a true solution which satisfies the conditions, and let the two-dimensional rectangular regionΩ=0,1×⌊0,1⌋; the true solution is as follows:(151)ux,y=u1x,y,u2x,y,u1x,y=10x2x−12yy−12y−1,u2x,y=−10xx−12x−1y2y−12,px,y=102x−12y−1.
It is proved that the velocity term of this vector function belongs to the nondivergence space, i.e.,∇⋅u=0, and the pressure term satisfies ∫Ωpxdx=0. By substituting ux,y,px,y into Navier–Stokes equation (1), the right-hand term of the equation can be obtained: fx,y=f1x,y,f2x,y. Since the viscosity coefficient has little effect on the accuracy of the equation, we only consider the case when the viscosity coefficient is ν=1. In this paper, we present the numerical results of the TNRFE method and the NRFE method when the value of h is different.
When the grid scale ish=1/25, the true solution and the calculation solution of the velocity component u1 are shown in Figure 12.
The contour of the true solution and the calculation solution of the velocity componentu1 is shown in Figure 13.
The true solution and the calculation solution of the velocity componentu2 are shown in Figure 14.
The contour of the true solution and the calculation solution of the velocity componentu2 is shown in Figure 15.
The vector of the true solution and the calculation solution of the velocity is shown in Figure16.
The images of the true solution and the calculation solution of the pressure are shown in Figure17.
The contours of the true solution and the calculation solution of the pressure are shown in Figure18.
From Tables2 and 3, we can see that the NRFE method and the TNRFE method have the same order convergence speed when h=OH2 for solving the Navier–Stokes equations. The TNRFE method we proposed is less time-consuming than the NRFE method; actually, it nearly saves half the CPU time. We plot the velocity by tracing diagonal points. When the grid scale is H=1/5andh=1/25, the true solution and the computational solution of velocity component u1 are represented in Figure 19.
The contour of the true solution and the calculation solution of the velocity componentu1 is shown in Figure 20.
The true solution and the calculation solution of the velocity componentu2 are shown in Figure 21.
The contour of the true solution and the calculated solution of the velocity componentu2 is shown in Figure 22.
The vector of the true solution and the calculated solution of the velocity is shown in Figure23.
The images of the true solution and the calculated solution of the pressure are shown in Figure24.
The contours of the true solution and the calculated solution of the pressure are shown in Figure25.Figure 12
True solution and calculation solution of velocity componentu1 by the NRFE method for the Navier–Stokes equations. (a) True solution of velocity component u1. (b) Calculation solution of velocity component u1.
(a)(b)Figure 13
True solution contour and calculation solution contour of velocity componentu1 by the NRFE method for the Navier–Stokes equations. (a) True solution contour of velocity component u1. (b) Calculation solution contour of velocity component u1.
(a)(b)Figure 14
True solution and calculation solution of velocity componentu2 by the NRFE method for the Navier–Stokes equations. (a) True solution of velocity component u2. (b) Calculation solution of velocity component u2.
(a)(b)Figure 15
True solution contour and calculation solution contour of velocity componentu2 by the NRFE method for the Navier–Stokes equations. (a) True solution contour of velocity component u2. (b) Calculation solution contour of velocity component u2.
(a)(b)Figure 16
True solution vector and calculation solution vector of velocity by the NRFE method for the Navier–Stokes equations. (a) True solution vector of velocity. (b) Calculation solution vector of velocity.
(a)(b)Figure 17
True solution and calculation solution of pressure by the NRFE method for the Navier–Stokes equations. (a) True solution of pressure. (b) Calculation solution of pressure.
(a)(b)Figure 18
True solution contour and calculation solution contour of pressure by the NRFE method for the Navier–Stokes equations. (a) True solution contour of pressure. (b) Calculation solution contour of pressure.
(a)(b)Table 2
Convergence study of the NRFE method for the Navier–Stokes equations in Example2.
hCPU (s)‖u−uh‖1,h/‖u‖1,h‖u−uh‖0/‖u‖0‖p−ph‖0/‖p‖01/910.052852646680.005120709900.006606373771/16290.050628430520.004692047150.006308510301/253910.048581969770.004314868840.00603642962Table 3
Convergence study of the TNRFE method for the Navier–Stokes equations in Example2.
HhCPU (s)‖u−uh‖1,h/‖u‖1,h‖u−uh‖0/‖u‖0‖p−ph‖0/‖p‖01/31/910.044944063710.006826184930.007697328201/41/16160.043320360000.006173997370.007295627301/51/251960.041808970720.005610628780.00693386430Figure 19
True solution and calculation solution of velocity componentu1 by the TNRFE method for the Navier–Stokes equations. (a) True solution of velocity component u1. (b) Calculation solution of velocity component u1.
(a)(b)Figure 20
True solution contour and calculation solution contour of velocity componentu1 by the TNRFE method for the Navier–Stokes equations. (a) True solution contour of velocity component u1. (b) Calculation solution contour of velocity component u1.
(a)(b)Figure 21
True solution and calculation solution of velocity componentu2 by the TNRFE method for the Navier–Stokes equations. (a) True solution of velocity component u2. (b) Calculation solution of velocity component u2.
(a)(b)Figure 22
True solution contour and calculation solution contour of velocity componentu2 by the TNRFE method for the Navier–Stokes equations. (a) True solution contour of velocity component u2. (b) Calculation solution contour of velocity component u2.
(a)(b)Figure 23
True solution vector and calculation solution vector of velocity by the TNRFE method for the Navier–Stokes equations. (a) True solution vector of velocity. (b) Calculation solution vector of velocity.
(a)(b)Figure 24
True solution and calculation solution of pressure by the TNRFE method for the Navier–Stokes equations. (a) True solution of pressure. (b) Calculation solution of pressure.
(a)(b)Figure 25
True solution contour and calculation solution contour of pressure by the TNRFE method for the Navier–Stokes equations. (a) True solution contour of pressure. (b) Calculation solution contour of pressure.
(a)(b)


## 7. Conclusion

The main work of this paper is as follows:(1)
The inequalities on the elements used to prove the standard polynomial approximation results are given.(2)
The existence, uniqueness, and convergence rate of the NRFE method for the steady Stokes equations and the Navier–Stokes equations are presented. The program design and numerical experiment results are given to verify the theory that the method has the same order of convergence speed as the conforming finite element (CFE) method.(3)
The stability and convergence rate of the Newton TNRFE method for the steady Navier–Stokes equations are proved, and the numerical results of the method are given.Further work is as follows:(1)
The NRFE method and the TNRFE method are extended to the unsteady Navier–Stokes equations, and corresponding numerical results are obtained.(2)
The rapid development of engineering and materials science in recent years has put forward new and higher requirements for the reliability of numerical models (i.e., the existence and uniqueness of solutions, convergence, and adaptability to the change of computational background) and nonlinear behavior. The problem of overall optimization of the numerical performance of the element is presented.



---

*Source: 7849163-2022-02-23.xml*