# On the Error Estimates of a New Operator  Splitting Scheme for the Navier-Stokes Equations with Coriolis Force

Authors: Chunjiang Shuai; Kaimin Teng; Hongen Jia
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

An operator splitting scheme is introduced for the numerical solution of the incompressible Navier-Stokes equations with Coriolis force. Under some mild regularity assumptions
on the continuous solution, error estimates and the stability analysis for the velocity and the pressure of the new operator splitting scheme are obtained. Some numerical results are presented to verify the
theoretical predictions.


---

## Body

## 1. Introduction

In this paper, we consider the numerical approximation of the unsteady Navier-Stokes equations with Coriolis force:(1.1)∂u∂t+(u·∇)u+2ω×u-νΔu+∇p=f,∇·u=0,
where Ω is a bounded region in Rd (d=2,3) with a sufficiently regular boundary ∂Ω. u is the velocity field, p is the pressure divided by the density (i.e., the kinematic pressure), ν=1/Re is the kinematic viscosity coefficients, Re is the Reynolds number, f is the vector of body forces, ω is the angular velocity vector, r is the radius vector from the center of coordinates, and 2ω×u is the so-called Coriolis force.For the sake of completeness, the equations should be supplemented with appropriate initial and boundary condition:(1.2)u(x,0)=u0(x)x∈Ω,u(x,t)=g(t,x)x∈∂Ω,t∈[0,T].
The difficulties for the numerical simulation of the incompressible flows are mainly of two kinds: nonlinearity and incompressibility. The velocity and the pressure are coupled by the incompressibility constraint, which requires that the solution spaces, to which the velocity and the pressure belong, verify the so-called inf-sup condition. To overcome these difficulties, operator splitting method, which can be viewed as the fractional step method, is introduced. Fractional step methods allow to separate the effects of the different operators appearing in the equation by splitting the time advancement into a series of substeps. In addition, the cost of simulation can be also reduced by using the fractional step method. However, these methods have a main disadvantage that splitting error is inevitable unless the operator is commute.The origin of this category of methods is contributed to the work of Chorin and Témam [1, 2]. They developed the so-called projection method, in which the second step consists of the projection of an intermediate velocity field onto the space of solenoidal vector field. The most attractive feature of projection methods is that, at each time step, one only needs to solve a sequence of decoupled elliptic equations for the velocity and the pressure, which makes it very efficient for large-scale numerical simulations. However, several issues related to these methods still deserve further analysis, and perhaps the most salient of these are the behavior of the computed pressure near boundaries and the stability of the pressure itself. The incompatibility of the projection boundary conditions may introduce a numerical boundary layer of size O(νΔt) [3, 4], where ν is the kinematic viscosity and Δt is the time step size. The end-of-step velocities of the projection do not converge in the space H01(Ω)d, since they do not satisfy the correct boundary conditions.These methods have been widely investigated.  Guermond et al. in [5] review theoretical and numerical convergence results available for projection methods. In [6–9], the analysis on first-order accurate schemes in the time size is presented. In [10, 11], Shen derived a second-order error estimates for the projection method.  Olshanskii et al. [12] proposed a projection method for the Navier-Stokes equations with Coriolis force and study the accuracy of its semidiscrete form. In [13], a new discrete projection method for the numerical solution of the Navier-Stokes equations with Coriolis force is presented, where the scheme is treated as an incomplete LU factorization of the transition operator for fully implicit time discretization. In [14], complex 3D simulations of the Stirred Tank Reactor model by a modified discrete projection method for the rotating incompressible flow are presented. Numerical experiments from [13, 14] show that including ω-term in the second step enhances stability and accuracy of the scheme for the case of dominating Coriolis forces.In this paper, we will consider the unsteady Navier-Stokes equations with Coriolis force (1.1). Using the technique developed in [7, 11] for the case of ω=0, a new scheme is introduced, which is a two-step scheme and allows to enforce the original boundary conditions of the problem in all substeps of the scheme. Some error estimates of both velocity and pressure for the proposed operator splitting scheme are given, which leads to the convergence of both the intermediate and the end-of-step velocities of the method to a continuous solution in the spaces L2(Ω)d and H01(Ω)d as in [15].The remainder of this paper is organized as follows. In Section2, we introduce some notations and assumptions, such as the regularity assumption for their solution. In Section 3, we present a new operator splitting scheme. In Section 4, the stability analysis is presented, then in Section 5, some error estimates for the intermediate, end-of-step velocity, and the pressure are given. Finally, in Section 6, some numerical results are presented to illustrate the theoretical results.


## 2. Function Setting

In order to study approximation scheme for problem (1.1), the following notations and assumptions are presented. we denote by (·,·) and ∥·∥ the inner product and norm in L2(Ω) or L2(Ω)d. The spaces H01(Ω) and H01(Ω)d are equipped with their usual norm; that is,
(2.1)∥u∥12=∫Ω‍|∇u(x)|2dx.
The norm in Hs(Ω) will be denoted simply by ∥·∥s. We will use 〈·,·〉 to denote the duality between H-s(Ω) and H0s(Ω) for all s>0.The following subspace is also introduced:(2.2)V={u∈H01(Ω)d:divu=0},H={u∈L2(Ω)d:divu=0,u·n=0}.
For the treatment of the convective term, the following trilinear form is given:
(2.3)b(u,v,w)=∫Ω‍(u·∇)v·wdx.
It is well known that b(·,·,·) is continuous in Hm1(Ω)×Hm2+1(Ω)×Hm3(Ω), provided m1+m2+m3≥d/2 if mi≠d/2,i=1,2,3, and this form is skew-symmetric with respect to its last two arguments, that is,
(2.4)b(u,v,w)=-b(u,w,v)∀u∈H,v,w∈H01(Ω)d.
In particular, we have
(2.5)b(u,v,v)=0∀u∈H,v∈H01(Ω)d,
and for d≤4,
(2.6)b(u,v,w)≤{c∥u∥1∥v∥1∥w∥1,c∥u∥∥v∥2∥w∥1,c∥u∥1∥v∥2∥w∥,c∥u∥∥v∥1∥w∥2,c∥u∥2∥v∥1∥w∥,c∥u∥1∥v∥1∥w∥1/2∥w∥11/2,c∥u∥1/2∥u∥11/2∥v∥1∥w∥1.
We also define the Stokes operator:
(2.7)Au=-PHΔu,∀u∈D(A)=V∩H2(Ω)d,
where PH is an orthogonal projector in the Hilbert space L2(Ω)d onto its subspace H. The Stokes operator A is an unbounded positive self-adjoint closed operator in H with domain D(A), and its inverse A-1 is compact in H. Having the following properties: there exists constant c1,c2>0, such that ∀u∈H,
(2.8)∥A-1u∥s≤c1∥u∥s-2fors=1,2,c2∥u∥-12≤(A-1u,u)≤c12∥u∥-12.
Furthermore, from (2.8), we will use (A-1u,u)1/2 as an equivalent norm of H-1(Ω)d for u∈H.For the purpose of this paper, we also need the following regularity assumptions:(A1)u0∈H1(Ω)d∩V,f∈L∞(0,T;L2(Ω)d)∩L2(0,T;H1(Ω)d). In the three-dimension case, we assume additionally (A2)supt∈[0,T]∥u(t)∥1≤M, where (A2) is automatically satisfied with some appropriate constant M when d=2.Under the regularity assumption (A1)-(A2), one can show that [9] (a)supt∈[0,T]{∥u(t)∥2+∥ut(t)∥+∥∇p(t)∥}≤M1,(b)∫0T‍∥ut(t)∥12≤M1.In addition, if we also assume that(A3)ft∈L2(0,T;L2(Ω))d holds, we have (c)∫0T‍∥utt∥-12dt≤M1which will be used in the sequel. Indeed, the estimates (a-b) and the estimate (c) were proved for the Navier-Stokes without Coriolis term in [9, 16], respectively. However, adding linear skew-symmetric term ω×u to the momentum equation does not change arguments from [16], but leads to (a)–(c) with constant M1 depending, in general, on Ω   [12]. Next, we cite the following lemma, which will be frequently used.Lemma 2.12.1 (Discrete Gronwall Lemma).
Letyn, hn, gn, and fn be nonnegative sequences satisfying
(2.9)ym+Δt∑n=0n=m‍hn≤B+Δt∑n=0n=m‍(gnyn+fn),withΔt∑n=0n=[T/Δt]‍gn≤M,∀0≤m≤[TΔt].
Assume Δtgn<1 and let σ=max0≤n≤[T/Δt](1-Δtn)-1, then
(2.10)ym+Δt∑n=0n=m‍hn≤exp(σM)(B+Δt∑n=0n=m‍fn)∀m≤[TΔt].
Hereafter, we will use c to denote a generic constant which depends only on Ω,ν,T, and constants from various Sobolev inequalities. We will denote M as a generic positive constant which may additionally depend on u0, f, ω.


## 3. New Operator Splitting Scheme

Equation (1.1) can be written as
(3.1)∂u∂t=A1+A2,
such that
(3.2)A1=-(u·∇)u+12νΔu,A2=12νΔu-∇p-2ω×u+f.
So an algorithm can be formulated as follows: for t∈[tn,tn+1](3.3){∂u~∂t=A1,u~(tn,x)=u(tn,x),→{∂u^∂t=A2,u^(tn,x)=u~(tn+1,x),
and take u(tn+1,x)≈u^ as the approximate solution of (1.1) at time tn+1.The scheme (3.3) has an irreducible splitting error of order O(Δt). Hence, using a higher-order time stepping scheme does not improve the overall accuracy. So a first-order accurate semidiscrete version can be obtained as follows: let u0=u0, we solve successively u~n+1 and {u^n+1,p^n+1} by
(3.4)u~n+1-u^nΔt+(u^n·∇)u~n+1-12νΔu~n+1=0,u~n+1|∂Ω=g(tn+1,x),(3.5)u^n+1-u~n+1Δt-12νΔu^n+1+∇p^n+1+2ω×u^n+1=f(tn+1),divu^n+1=0,u^n+1=g(tn+1,x).
Note that we have omitted the dependency to x of the function f to simplify our notations; we will do so for {u,p}.As can be seen in (3.5), the main difference between this method and the standard projection method is the introduction of a viscous term in the incompressible step, which allows the imposition of the original boundary condition (2.6) on the end-of-step velocity u^n+1. Similar ideas can be found in the θ-scheme [17] in which viscosity and incompressibility are also coupled, and some other methods such as [18–20], all of which involve an incompressible step with part of the viscous term. It leads to convergence of both the intermediate and the end-of-step velocities of the method to a continuous solution in space L2(Ω) and H1(Ω). In comparison with the θ scheme, our scheme is two steps instead of three steps in θ scheme. Moreover, the fact that un+1 satisfies the correct boundary conditions will allow us to obtain improved error estimates comparative with the standard projection method.Denoting the corresponding right-hand side byf, at each time step, we have to solve the following two subproblems:
(3.6)αu+(w·∇)u+βΔu=f,u~n+1|∂Ω=g(tn+1,x),(3.7)γu+βΔu+∇p=f,divu^n+1=0,u^n+1=g(tn+1,x),
with α=1/Δt,β=-(1/2)ν,γ=1/Δt+2×ω.The first step of the method is a linearized elliptic problem, which can be seen as a linear Burger’s problem. The second step is a generalized Stokes problem. To solve problem (3.6), the fixed point iterative technique is used as in [21, 22], which is cheaper than the conjugate gradient method used by the least-square technique in the corresponding advective subproblems appearing in Glowinski’s θ-scheme. To solve problem (3.7), the efficient technique of the functional equation satisfied by the pressure is used, which is used in the Glowinski’s θ-scheme for the corresponding Stokes problem; that is, conjugate gradient method is applied on the variational formulation of such an equation. One defect of this method is that the discrete inf-sup compatibility condition should be satisfied.


## 4. Stability Analysis

For the sake of simplicity, we will only consider the homogeneous boundary conditionu(x,t)|∂Ω=0, that is, g(t,x)=0 for the scheme (3.4)-(3.5).Theorem 4.1.
Under the assumptions(A1)-(A2), there exists a constant c, such that
(4.1)u~n∈L2(0,T,H1(Ω)d)∩L∞(0,T,L2(Ω)d),u^n∈L2(0,T,H1(Ω)d)∩L∞(0,T,L2(Ω)d).Proof.
We take the inner product of (3.4) with 2Δtu~n+1 to get
(4.2)∥u~n+1-u^n∥2+∥u~n+1∥2-∥u^n∥2+νΔt∥u~n+1∥12=0.
Next, taking the inner product of (3.5) with 2Δtu^n+1, and using the condition divu^n+1=0, we obtain
(4.3)∥u~n+1-u^n+1∥2+∥u^n+1∥2-∥u~n+1∥2+νΔt∥u^n+1∥12=2Δt(f(tn+1),u^n+1).
Combing (4.2) with (4.3), and using the Young’s inequality, we have
(4.4)∥u^n+1∥2+∥u~n+1-u^n+1∥2+∥u~n+1-u^n∥2+νΔt(∥u~n+1∥12+∥u~n+1∥12)≤cΔt∥(f(tn+1)∥2+Δt∥u^n+1∥2+∥u^n∥2.
Summing up the inequality (4.4) for n=0,…,r≤N, we obtain
(4.5)∥u^r+1∥2+∑n=0n=r‍∥u~n+1-u^n+1∥2+∑n=0n=r‍∥u~n+1-u^n∥2+Δtν∑n=0n=r‍∥u^n+1∥12+νΔt∑n=0n=r‍∥u~n+1∥12≤cΔt∑n=0n=r‍∥f(tn+1)∥02+∥u0∥2+Δt∑n=0n=r‍∥u^n+1∥2.
Applying the discrete Gronwall lemma to the above inequality, we obtain
(4.6)∥u^r+1∥2+∑n=0n=r‍∥u~n+1-u^n+1∥2+∑n=0n=r‍∥u~n+1-u^n∥2+Δtν∑n=0n=r‍∥u^n+1∥12+νΔt∑n=0n=r‍∥u~n+1∥12≤cΔt∑n=0n=r‍∥f(tn+1)∥02+∥u0∥2.
Thus, using the regularity properties of the continuous solution u, for arbitrary n, we have
(4.7)∥u^n+1∥2+∑n=0n=N‍(∥u~n+1-u^n+1∥2+∥u~n+1-u^n∥2)+νΔt∑n=0n=N‍(∥u^n+1∥12+∥u~n+1∥12)≤c,
which means that
(4.8)u~n∈L2(0,T,H1(Ω)d),u^n∈L2(0,T,H1(Ω)d)∩L∞(0,T,L2(Ω)d).
From (4.2), it yields(4.9)∥u~n+1∥2≤∥u^n∥2.
So we have u~n∈L∞(0,T,L2(Ω)d). The proof is completed.Remark 4.2.
The formula (4.7) can be viewed as the discrete version of the classical energy estimate for the Navier-Stokes equations [23]. From (4.6), we have
(4.10)∥u^N+1∥2≤cΔt∑n=0n=N‍∥f(tn+1)∥02+∥u0∥2.
This estimate provides a meaningful bound for ∥u^N+1∥2 for the first few time steps, that is, for small T.


## 5. Error Analysis

In this section, we present an error analysis of the operator splitting scheme introduced in the previous section. Firstly, we define the semidiscrete velocity error as(5.1)e^n+1=u(tn+1)-u^n+1,e~n+1=u(tn+1)-u~n+1.
We give a first estimate for e^n+1 and e~n+1 which shows that u^n+1 and u~n+1 are both order 1/2 approximations to u in l∞(L2(Ω)d) and l2(H1(Ω)d).Theorem 5.1.
Under the regularity assumptions(A1)–(A3), there exists a constant M, such that
(5.2)∥e^N+1∥2+∥e~N+1∥2+∑n=0n=N‍(∥e~n+1-e^n∥2+∥e^n+1-e~n+1∥2+νΔt(∥e^n+1∥12+∥e~n+1∥12))≤MΔt.Proof.
LetRn be the truncation error defined by
(5.3)1Δt(u(tn+1)-u(tn))-νΔu(tn+1)+2ω×u(tn+1)+(u(tn+1)·∇)u(tn+1)+∇p(tn+1)=f(tn+1)+Rn,
where Rn is the integral residual of the Taylor series, that is,
(5.4)Rn=1Δt∫tntn+1‍(t-tn)uttdt.
By subtracting (3.4) from (5.3), we obtain
(5.5)e~n+1-e^nΔt-12νΔe~n+1-12νΔu(tn+1)+2ω×u(tn+1)=(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)+Rn-∇p(tn+1)+f(tn+1).
The nonlinear terms on the right-side can be split into three terms:
(5.6)(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)=-(e^n·∇)u~n+1+((u(tn)-u(tn+1))·∇)u~n+1-(u(tn+1)·∇)e~n+1.
Taking the inner product of (5.5) with 2Δte~n+1, using the identity (a-b,2a) =|a|2+|a+b|2-|b|2, we obtain
(5.7)∥e~n+1∥2+∥e~n+1-e^n∥2-∥e^n∥2+νΔt∥e~n+1∥12-νΔt(Δu(tn+1),e~n+1)+4Δt(ω×u(tn+1),e~n+1)=2Δt<Rn,e~n+1>-2Δt(∇p(tn+1),e~n+1)-2Δtb(e^n,u~n+1,e~n+1)+2Δtb(u(tn)-u(tn+1),u~n+1,e~n+1)-2Δtb(u(tn+1),e~n+1,e~n+1)+2Δt(f(tn+1),e~n+1).
On the other hand, we derive from (3.5) that
(5.8)e^n+1-e~n+1Δt+12νΔu^n+1-∇p^n+1-2ω×u^n+1+f(tn+1)=0.
Taking the inner product of the last equality with 2Δte^n+1, and using dive^n+1=0, we obtain
(5.9)∥e^n+1∥2+∥e^n+1-e~n+1∥2-∥e~n+1∥2+Δtν(Δu^n+1,e^n+1)-4Δt(ω×u^n+1,e^n+1)+2Δt(f(tn+1),e^n+1)=0.
Combing (5.7) with (5.9), we obtain
(5.10)∥e^n+1∥2+∥e~n+1-e^n∥2+∥e^n+1-e~n+1∥2+Δtν∥e^n+1∥12+νΔt∥e~n+1∥12=2Δt<Rn,e~n+1>-2Δt(∇p(tn+1),e~n+1)-2Δtb(e^n,u~n+1,e~n+1)+2Δtb(u(tn)-u(tn+1),u~n+1,e~n+1)-2Δtb(u(tn+1),e~n+1,e~n+1)+2Δt(f(tn+1),e~n+1-e^n+1)+Δtν(Δu(tn+1),e~n+1-e^n+1)+∥e^n∥2-4Δt(ω×u(tn+1),e~n+1-e^n+1).
We bound each term in the right-hand side of (5.10) independently.
Taylor residual term:(5.11)2Δt〈Rn,e~n+1〉≤2Δt∥Rn∥-1∥e~n+1∥1≤Δtν6∥e~n+1∥12+cΔt2∫tntn+1‍∥utt∥-12dt.
Pressure gradient term:
(5.12)-2Δt(∇p(tn+1),e~n+1)=-2Δt(∇p(tn+1),e~n+1-e^n)≤12∥e~n+1-e^n∥2+cΔt2∥∇p(tn+1)∥2.
Nonlinear term:
(5.13)-2Δtb(e^n,u~n+1,e~n+1)=-2Δtb(e^n,u(tn+1),e~n+1)≤cΔt∥e^n∥∥e~n+1∥1∥u(tn+1)∥2≤Δtν6∥e~n+1∥12+cΔt∥e^n∥2,2Δtb(u(tn)-u(tn+1),u~n+1,e~n+1)=-2Δtb(u(tn)-u(tn+1),e~n+1,u~n+1)=-2Δtb(u(tn)-u(tn+1),e~n+1,u(tn+1))≤cΔt∥u(tn)-u(tn+1)∥∥e~n+1∥1∥u(tn+1)∥2≤Δtν6∥e~n+1∥12+cΔt2∫tntn+1‍∥ut∥2dt,-2Δtb(u(tn+1),e~n+1,e~n+1)=0.
The external term:
(5.14)2Δt(f(tn+1),e~n+1-e^n+1)≤16∥e~n+1-e^n+1∥2+cΔt2∥f(tn+1)∥2.
The viscous term:
(5.15)2Δt(νΔu(tn+1),e~n+1-e^n+1)≤16∥e~n+1-e^n+1∥2+cΔt2∥Δu(tn+1)∥2≤16∥e~n+1-e^n+1∥2+cΔt2∥u(tn+1)∥22.
The rotating term:
(5.16)-4Δt(ω×u(tn+1),e~n+1-e^n+1)≤cΔt2∥ω×u(tn+1)∥2+16∥e~n+1-e^n+1∥2≤cΔt2∥ω∥2∥u(tn+1)∥2+16∥e~n+1-e^n+1∥2.
Inserting the above estimates into (5.10), we obtain
(5.17)∥e^n+1∥2+12∥e~n+1-e^n∥2+12∥e^n+1-e~n+1∥2+νΔt∥e^n+1∥12+12νΔt∥e~n+1∥12≤cΔt2∫tntn+1‍∥utt∥-12dt+cΔt2∥∇p(tn+1)∥2+cΔt∥e^n∥2+cΔt2∫tntn+1‍∥ut∥2dt+cΔt2∥f(tn+1)∥2+cΔt2∥u(tn+1)∥22+cΔt2∥ω∥2∥u(tn+1)∥2.
Summing up the inequality (5.17) for n=0,…,N, we get
(5.18)∥e^N+1∥2+∑n=0n=N‍(12∥e~n+1-e^n∥2+12∥e^n+1-e~n+1∥2+νΔt∥e^n+1∥12+12νΔt∥e~n+1∥12)≤cΔt(Δt∫0tT‍∥utt∥-12dt+supt∈[0,T]∥∇p(tn+1)∥2+Δt∫0T‍∥ut∥2dt+∑n=0n=N‍Δt∥f(tn+1)∥2+∑n=0n=N‍Δt∥u(tn+1)∥2+∑n=0n=N‍Δt∥ω∥2∥u(tn+1)∥2).
By applying the discrete Gronwall lemma to the above inequality, we derive
(5.19)∥e^N+1∥2+∑n=0n=N‍(12∥e~n+1-e^n∥2+12∥e^n+1-e~n+1∥2+νΔt∥e^n+1∥12+12νΔt∥e~n+1∥12)≤cΔt(Δt∫0T‍∥utt∥-12dt+supt∈[0,T]∥∇p(tn+1)∥2+Δt∫0T‍∥ut∥2dt+∑n=0n=N‍Δt∥f(tn+1)∥2+∑n=0n=N‍Δt∥u(tn+1)∥22+∑n=0n=N‍Δt∥ω∥2∥u(tn+1)∥2).
Using the regularity properties of u, we obtain
(5.20)∥e^N+1∥2+∑n=0n=N‍(∥e~n+1-e^n∥2+∥e^n+1-e~n+1∥2+νΔt∥e^n+1∥12+νΔt∥e~n+1∥12)≤MΔt.
Finally, the bounds for u~n+1 follow from (5.20) and the triangle inequality. Theorem 5.1  is proved.Remark 5.2.
Theorem5.1 shows that
(5.21)∥u~n+1∥1≤M,∥u^n+1∥1≤M,
since ∥e^n+1∥1≤M,∥e~n+1∥1≤M. Moreover, we also have
(5.22)∥u~n+1∥0≤MΔt1/2,∥u^n+1∥0≤MΔt1/2.
Next, we will use the previous result to improve the error estimates for the velocity and give an error estimate for the pressure as well. The result shows that bothu^n+1 and u~n+1 are weakly first-order approximations to u(tn+1) in l2(L2(Ω)d).Theorem 5.3.
Under the regularity assumptions(A1)–(A3), there exists a constant M, such that
(5.23)∥e^N+1∥-12+∑N=0N‍(∥e^n+1-e^n∥-12+νΔt∥e~n+1∥2+νΔt∥e^n+1∥2)≤MΔt2,(5.24)Δt∑n=0n=N‍∥p(tn+1)-p^n+1∥L02Ω≤MΔt.Proof.
Taking the sum of (3.4) and (3.5), we obtain
(5.25)u^n+1-u^nΔt+(u^n·∇)u~n+1-12νΔu~n+1-12νΔu^n+1+∇p^n+1+2ω×u^n+1=f(tn+1).
Let us denote
(5.26)q^n+1=p(tn+1)-p^n+1.
Subtracting (5.25) from (5.3), we obtain
(5.27)e^n+1-e^nΔt-12νΔe~n+1-12νΔe^n+1+∇q^n+1+2ω×e^n+1=(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)+Rn.Taking the inner product of the expression (5.27) with 2ΔtA-1e^n+1, we obtain
(5.28)(e^n+1,A-1e^n+1)-(e^n,A-1e^n)+(e^n+1-e^n,A-1(e^n+1-e^n))-νΔt(Δe~n+1,A-1e^n+1)-νΔt(Δe^n+1,A-1e^n+1)+4Δt(ω×e^n+1,A-1e^n+1)=2Δtb(u^n,u~n+1,A-1e^n+1)-2Δtb(u(tn+1),u(tn+1),A-1e^n+1)+2Δt(Rn,A-1e^n+1).
Since
(5.29)-νΔt(Δe~n+1,A-1e^n+1)≥νΔt(14∥e^n+1∥2-c∥e~n+1-e^n+1∥2)-νΔt(Δe^n+1,A-1e^n+1)=νΔt∥e^n+1∥2,
together with (5.28), it yields
(5.30)∥e^n+1∥-12-∥e^n∥-12+∥e^n+1-e^n∥-12+54νΔt∥e^n+1∥2≤2Δtb(u^n,u~n+1,A-1e^n+1)-2Δtb(u(tn+1),u(tn+1),A-1e^n+1)+2Δt(Rn,A-1e^n+1)+cΔt∥e~n+1-e^n+1∥2-4Δt(ω×e^n+1,A-1e^n+1).
Similar to (5.6) for the nonlinear term, together with (5.30), it yields
(5.31)∥e^n+1∥-12-∥e^n∥-12+∥e^n+1-e^n∥-12+54νΔt∥e^n+1∥2=-2Δtb(e^n,u~n+1,A-1e^n+1)+2Δtb(u(tn)-u(tn+1),u~n+1,A-1e^n+1)-2Δtb(u(tn+1),e~n+1,A-1e^n+1)+2Δt(Rn,A-1e^n+1)+cΔt∥e~n+1-e^n+1∥2-4Δt(ω×e^n+1,A-1e^n+1).
We will focus on the right-hand side as follows.
The Coriolis term is estimated as follows:(5.32)∥4Δt(ω×e^n+1,A-1e^n+1)∥≤cΔt∥A-1e^n+1∥1∥e^n+1∥≤cΔt∥e^n+1∥-1∥e^n+1∥≤νΔt20∥e^n+1∥2+cΔt∥e^n+1∥-12.
For the Taylor residual term, we have
(5.33)2Δt(Rn,A-1e^n+1)≤cΔt∥Rn∥-1∥A-1e^n+1∥1≤Δt∥e^n+1∥-12+cΔt2∫tntn+1‍∥utt∥-12dt.
For the nonlinear term, it yields
(5.34)-2Δtb(e^n,u~n+1,A-1e^n+1)=2Δtb(e^n,A-1e^n+1,u~n+1)=2Δtb(e^n,A-1e^n+1,u(tn+1))-2Δtb(e^n,A-1e^n+1,e~n+1)=T1+T2.
Using (2.6)-(2.8), we derive
(5.35)T1≤cΔt∥e^n∥∥A-1e^n+1∥1∥u(tn+1)∥2≤cΔt∥e^n∥∥e^n+1∥-1≤cΔt(∥e^n+1∥+∥e^n+1-e~n+1∥+∥e~n+1-e^n∥)∥e^n+1∥-1≤νΔt20∥e^n+1∥2+cΔt(∥e^n+1-e~n+1∥2+∥e~n+1-e^n∥2)+cΔt∥e^n+1∥-12,
and together with Theorem 5.1,
(5.36)T2≤cΔt∥e^n∥∥A-1e^n+1∥2∥e~n+1∥1≤cΔt∥e^n∥∥e^n+1∥∥e~n+1∥1≤cΔt3/2∥e^n+1∥∥e~n+1∥1≤νΔt20∥e^n+1∥2+cΔt2∥e~n+1∥1,-2Δtb(u(tn+1),e~n+1,A-1e^n+1)=2Δtb(u(tn+1),A-1e^n+1,e~n+1)≤cΔt∥u(tn+1)∥2∥A-1e^n+1∥1∥e~n+1∥≤cΔt∥e^n+1∥-1∥e~n+1∥≤cΔt∥e^n+1∥-1(∥e^n+1-e~n+1∥+∥e^n+1∥)≤νΔt20∥e^n+1∥2+cΔt∥e^n+1-e~n+1∥2+cΔt∥e^n+1∥-12.
Similarly,
(5.37)2Δtb(u(tn)-u(tn+1),u~n+1,A-1e^n+1)≤∥u(tn)-u(tn+1)∥∥u~n+1∥1∥A-1e^n+1∥2≤νΔt16∥e^n+1∥2+cΔt2∫tntn+1‍∥ut∥2dt.
Inserting the above inequality into (5.31), we obtain
(5.38)∥e^n+1∥-12-∥e^n∥-12+∥e^n+1-e^n∥-12+νΔt∥e^n+1∥2=cΔt2(∫tntn+1‍∥utt∥-12dt+∫tntn+1‍∥ut∥2dt)+cΔt(∥e^n+1-e~n+1∥2+∥e~n+1-e^n∥2)+cΔt2∥e~n+1∥1+cΔt∥e^n+1∥-12.
Taking the sum of the above inequality for n=0,…,N, using the regularity assumption of the solution u and Theorem 5.1, yields
(5.39)∥e^N+1∥-12+∑N=0N‍(∥e^n+1-e^n∥-12+νΔt∥e^n+1∥2)≤mΔt2+cΔt∑N=0N‍∥e^n+1∥-12.
By applying the discrete Gronwall lemma to the last inequality, we obtain
(5.40)∥e^N+1∥-12+∑N=0N‍(∥e^n+1-e^n∥-12+νΔt∥e^n+1∥2)≤MΔt2.
For u~n+1, we have
(5.41)νΔt∑N=0N‍∥e~n+1∥2≤νΔt∑N=0N‍(∥e^n+1∥+∥e~n+1-e^n+1∥)≤MΔt2,
and together with (5.40), we derive (5.23).
Next, we derive the estimate for the pressure; we recast (5.27) as
(5.42)∇q^n+1=12νΔe~n+1+12νΔe^n+1-e^n+1-e^nΔt+(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)+Rn-2ω×e^n+1,
firstly, by using (2.6) and Theorem 5.3, for all v∈H01(Ω)d,
(5.43)(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)≤c∥u(tn+1)-u(tn)∥∥u(tn+1)∥2∥v∥1+c∥e^n∥1∥u(tn+1)∥1∥v∥1+c∥u^n∥1∥e~n+1∥1∥v∥1≤(∥e~n+1∥1+∥e^n∥1+∥u(tn+1)-u(tn)∥)∥v∥1.
Using the Schwarz inequality, we have also, for all v∈H01(Ω)d,
(5.44)(12νΔe~n+1+12νΔe^n+1-e^n+1-e^nΔt-2ω×e^n+1+Rn,v)≤(1Δt∥e^n+1-e^n∥-1+∥Rn∥-1+12ν∥e^n+1∥1+12ν∥e~n+1∥1+c∥e^n+1∥1)∥v∥1.
Finally, we derive
(5.45)n∥q^n+1∥L02Ω≤csupv∈H01(Ω)d(∇q^n+1,v)∥v∥1≤cΔt∥e^n+1-e^n∥-1+c(∥Rn∥-1+∥e~n+1∥1+∥e^n+1∥1+∥e^n∥1+∥u(tn+1)-u(tn)∥).
Therefore, by using Theorem 5.1 and (5.40), we obtain (5.24).The error estimate of Theorem5.3 can be improved to first order on the norms of l∞(L2(Ω)d) and l2(H01(Ω)d) for the end-of-step velocities u^n+1.Theorem 5.4.
Under the regularity assumptions(A1)–(A3), for small enough Δt, there exists a constant M, such that
(5.46)∥e^N+1∥2+∑n=0n=N‍∥e^n+1-e^n∥2+νΔt∑n=0n=N‍(∥e~n+1∥12+∥e^n+1∥12)≤MΔt2.Proof.
Taking the inner product of (5.27) with 2Δte^n+1, we have
(5.47)∥e^n+1∥2-∥e^n∥2+∥e^n+1-e^n∥2+νΔt∥e~n+1∥12+νΔt∥e^n+1∥12=2Δtb(u^n,u~n+1,e^n+1)-2Δtb(u(tn+1),u(tn+1),e^n+1)+2Δt〈Rn,e^n+1〉.
The estimates below are obtained on the right-hand term of (5.47):
(5.48)2Δt〈Rn,e^n+1〉≤νΔt10∥e^n+1∥12+cΔt2∫tntn+1‍∥utt∥-12dt.
For the nonlinear term, similar to (5.6), the below estimates are obtained:
(5.49)-2Δtb(u(tn+1),e~n+1,e^n+1)≤νΔt10∥e^n+1∥12+cΔt∥e~n+1∥2.
For the remainder of nonlinear term,  it yields
(5.50)2Δtb(u(tn)-u(tn+1),e~n+1,e^n+1)≤∥u(tn)-u(tn+1)∥1∥e~n+1∥1∥e^n+1∥1≤cΔt2∫tntn+1‍∥ut∥12dt+νΔt10∥e^n+1∥12,-2Δtb(e^n,u~n+1,e^n+1)=2Δtb(e^n,e~n+1,e^n+1)-2Δtb(e^n,u(tn+1),e^n+1).
Thus,
(5.51)2Δtb(e^n,e~n+1,e^n+1)≤cΔt∥e^n∥1∥e^n+1∥1∥e~n+1∥1/2∥e~n+1∥11/2≤cΔt∥e^n∥1∥e^n+1∥1∥e~n+1∥1/2≤cΔt5/4∥e^n∥1∥e^n+1∥1≤cΔt3/2∥e^n∥12+νΔt10∥e^n+1∥12,-2Δtb(e^n,u(tn+1),e^n+1)≤cΔt∥e^n∥0∥u(tn+1)∥2∥e^n+1∥1≤cΔt∥e^n∥0∥e^n+1∥1≤νΔt10∥e^n+1∥12+cΔt∥e^n∥2,
where we have used Theorem 5.1 and formula (2.6).
Taking the sum of the formula (5.47) for n from 0 to N, together with the above estimates, we get
(5.52)∥e^N+1∥2+∑n=0n=N‍∥e^n+1-e^n∥2+νΔt∑n=0n=N‍(∥e~n+1∥12+12∥e^n+1∥12)=CΔt2∫0T‍∥utt∥-12dt+cΔt∑n=0n=N‍∥e~n+1∥2+cΔt2∫0T‍∥ut∥12dt+cΔt3/2∑n=0n=N‍∥e^n∥12+cΔt∑n=0n=N‍∥e^n∥2.
By virtue of the formula (5.23) and the regularity assumption (A2), (A3), we obtain(5.53)∥e^N+1∥2+∑n=0n=N‍∥e^n+1-e^n∥2+νΔt∑n=0n=N‍(∥e~n+1∥12+12∥e^n+1∥12)=CΔt2+cΔt3/2∑n=0n=N‍∥e^n∥12+cΔt∑n=0n=N‍∥e^n+1∥2.
For sufficiently small Δt, we can take the last term to left side and apply the discrete Gronwall lemma to the last inequality, so the proof is completed.Theorem 5.5.
Under the regularity assumptions(A1)–(A3), for small enough Δt, there exists a positive constant δ, such that δ<1/8 and 4MΔt∥e^N∥2ν-1<1, then
(5.54)∑n=0N‍∥e~n+1-e~n∥2≤mΔt2.Proof.
We shift the indexn+1→n in (3.5) to get
(5.55)e^n-e~nΔt+12νΔu^n-∇p^n-2ω×u^n+f(tn)=0,
and taking the sum with (5.5), we obtain
(5.56)e~n+1-e~nΔt-ν2Δe~n+1=(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)+Rn-∇(p(tn+1)-p^n)+(f(tn+1)-f(tn))+ν2Δe^n+ν2Δ(u(tn+1)-u(tn))-2ω×(u(tn+1)-u(tn))-2ω×e^n.
Taking the inner product of (5.56) with Δt(e~n+1-e~n), the left-hand term of (5.56) can be written as
(5.57)∥e~n+1-e~n∥2+νΔt4(∥e~n+1∥12+∥e~n+1-e~n∥12-∥e~n∥12).
Now, we give the estimates of the right-hand term of (5.56):
(5.58)-2Δt(ω×(u(tn+1)-u(tn)),e~n+1-e~n)≤cΔt3∫tntn+1‍∥ut∥2dt+δ∥e~n+1-e~n∥2,-2Δt(ω×e^n,e~n+1-e~n)≤cΔt2∥e^n∥2+δ∥e~n+1-e~n∥2≤MΔt3+δ∥e~n+1-e~n∥2,
where we have used Theorem 5.1. Simultaneously,
(5.59)νΔt2(Δ(u(tn+1)-u(tn)),e~n+1-e~n)=-νΔt2(∇(u(tn+1)-u(tn)),∇(e~n+1-e~n))≤cΔt∥u(tn+1)-u(tn)∥1∥e~n+1-e~n∥1≤cΔt∥u(tn+1)-u(tn)∥12+νΔt8∥e~n+1-e~n∥12≤cΔt2∫tntn+1‍∥ut∥12dt+νΔt8∥e~n+1-e~n∥12,νΔt2(Δe^n,e~n+1-e~n)=-νΔt2(∇e^n,∇(e~n+1-e~n))≤cΔt(∥e^n∥1+∥e~n+1-e~n∥1)≤cΔt(∥e^n∥1+∥e~n+1∥1+∥e~n∥1)≤MΔt2,
where we have used Theorem 5.4, and
(5.60)Δt(f(tn+1)-f(tn),e~n+1-e~n)≤cΔt3∫tntn+1‍∥ft∥2dt+δ∥e~n+1-e~n∥2.
For the Taylor residual term, we have
(5.61)Δt(Rn,e~n+1-e~n)≤δ∥e~n+1-e~n∥2+cΔt2∫tntn+1‍t∥utt∥2dt.
For the pressure term, since dive^n+1=0, respectively dive^n=0, we obtain
(5.62)-Δt(∇(p(tn+1)-p^n),e~n+1-e~n)=Δt(p(tn+1)-p^n,∇(e~n-e~n+1))≤Δt∥p(tn+1)-p^n∥∥e~n+1-e~n∥1≤(∥p(tn+1)∥+∥p^n∥)(Δt∥e~n+1∥1+Δt∥e~n∥1)≤MΔt2,
where Theorems 5.3 and 5.4 are used.
For the trilinear term, we consider the below splitting:(5.63)(u^n·∇)u~n+1-(u(tn+1)·∇)u(tn+1)=e^n·∇e~n+1-(u(tn+1)-u(tn))·∇u~n+1-u(tn+1)·∇e~n+1-e^n·∇u(tn+1).
So, based on the formula (2.6) and regularity assumption, it yields
(5.64)Δtb(u(tn+1),e~n+1,e~n+1-e~n)≤δ∥e~n+1-e~n∥2+cΔt2∥e~n+1∥12,Δtb(u(tn+1)-u(tn),e~n+1,e~n+1-e~n)≤δ∥e~n+1-e~n∥2+cΔt3∥u~n+1∥1∫tntn+1‍∥ut∥22dt≤δ∥e~n+1-e~n∥2+MΔt3,Δtb(e^n,u(tn+1),e~n+1-e~n)≤δ∥e~n+1-e~n∥2+cΔt2∥e^n∥12,Δtb(e^n,e~n+1,e~n+1-e~n)≤cΔt∥e^n∥1∥e~n+1∥1∥e~n+1-e~n∥11/2∥e~n+1-e~n∥1/2≤cΔt3/2∥e^n∥12∥e~n+1∥12+Δtνδ∥e~n+1-e~n∥1∥e~n+1-e~n∥≤cΔt3/2∥e^n∥12∥e~n+1∥12+νΔt8∥e~n+1-e~n∥12+δ∥e~n+1-e~n∥2.
Using (5.56)–(5.64), the regularity assumption, and Theorem 5.1, we obtain
(5.65)(1-8δ)∥e~n+1-e~n∥2+νΔt4(∥e~n+1∥12-∥e~n∥12)≤m(Δt2+Δt3/2∥e^n∥12∥e~n+1∥12).
Summing up the above inequality for n=0,…,N, we have
(5.66)∑n=0N‍c0∥e~n+1-e~n∥2+νΔt4∥e~N+1∥12≤m(Δt2+∑n=0N‍Δt3/2∥e^n∥12∥e~n+1∥12),
where we assume 1-8δ>c0, c0 is a positive constant.
Now, we assume thatΔt is sufficiently small such that 4MΔt∥e^N∥2ν-1<1 holds. (Note that e^N is uniformly bounded due to Theorem 5.1, then by virtue of the discrete Gronwall inequality, the proof is completed).


## 6. Numerical Results

In this section, we will give some numerical results to verify the theoretical analysis for the new operator splitting method. We solve the system of the incompressible Navier-Stokes equations with Coriolis force term (1.1) with homogeneous Dirichlet boundary conditions on the velocity. The exact solution (u,p) is chosen as follows:
(6.1)u1=πsin(t)sin(2πy)sin2(πx),u2=-πsin(t)sin(2πx)sin2(πy),p=sin(t)cos(πx)sin(πy).
The initial condition is set equal to the exact solution and f is computed by evaluating the momentum equation of problem (1.1) for the exact solution.The following setting is chosen:Re=100, T=1, Δt={0.1,0.05,0.01,0.005}, |ω|=10 and the uniform mesh with the mesh-size 1/30. The stream line and pressure contours are shown by Figures 1 and 2. The experimental rates of convergence with respect to the time size Δt are given by Table 1, where
(6.2)Eu0=∥uanalyt(t)-unumer(t)∥L2(Ω)d,Eu1=(Δt∑n=1n=N‍∥uanalyt(nΔt)-unumer(nΔt)∥H1(Ω)d2)1/2,N=T/Δt, uanalyt(t) and panalyt(t) are gotten from (6.1), and unumer(t) and pnumer(t) are corresponding numerical values. From the graphics and the table above, one can observe that the numerical results are in good agreement with the theoretical analysis.Table 1
p2-p1 element with different time step.
Δ
t
Numerical solution
Convergence rate
E
u
0
E
u
1
E
u
0 rate
E
u
1 rate
0.100
0.18258700
0.41101800
/
/
0.050
0.09230720
0.14637600
0.98407
1.48950
0.010
0.01632060
0.01402730
1.07660
1.45710
0.005
0.00806472
0.00705089
1.01700
0.99236Exact solution:T=1.
(a)
(b)Numerical solution:T=1, dt=0.01.
(a)
(b)



---

*Source: 105735-2012-12-09.xml*