G. Zhou, Y. Zhou, B. Ahmad, a. Alsaedi: Finite difference/element method …

# Finite difference/element method for time-fractional Navier-Stokes equations

Guang-an Zou, Yong Zhou, Bashir Ahmad, Ahmed Alsaedi
###### Abstract.

We apply a composite idea of semi-discrete finite difference approximation in time and Galerkin finite element method in space to solve the Navier-Stokes equations with Caputo derivative of order . The stability properties and convergence error estimates for both the semi-discrete and fully discrete schemes are obtained. Numerical example is provided to illustrate the validity of theoretical results.

Keywords: Time-fractional Navier-Stokes equations; finite difference approximation; finite element method; error estimates; numerical examples

AMS Subject Classification: 76D05, 65N30, 65N12

School of Mathematics and Statistics, Henan University, Kaifeng 475004, P. R. China
Faculty of Mathematics and Computational Science, Xiangtan University
Xiangtan, Hunan 411105, P.R. China
Nonlinear Analysis and Applied Mathematics (NAAM) Research Group, Faculty of Science
King Abdulaziz University, Jeddah 21589, Saudi Arabia

## 1. Introduction

In this paper, we study the following Navier-Stokes equations with time-fractional derivative in a bounded subset of with a smooth boundary :

 ⎧⎪⎨⎪⎩CDαtu+(u⋅∇)u−νΔu+∇p=f,   div u=0, ∀ (x,t)∈Ω×(0,T], u(x,t)|∂Ω=0, t∈(0,T], u(x,0)=u0, x∈Ω, (1.1)

where represents the Caputo-type fractional derivative of order , denotes the velocity field at a point and , is viscosity coefficient, represents the pressure field, is the external force and is the initial velocity.

Notice that the problem (1.1) reduces to the classical Navier-Stokes equations (NSEs) for . The existence and non-existence of solutions for the NSEs have been discussed in [1]. Chemin et al.[2] studied the global regularity for the large solutions to the NSEs. Miura [3] focused on the uniqueness of mild solutions to the NSEs. Germain [4] presented the uniqueness criteria for the solutions of the Cauchy problem associated to the NSEs. The existence of global weak solutions for supercritical NSEs was discussed [5]. The lower bounds on blow up solutions for the NSEs in homogeneous Sobolev spaces were studied in [6]. The numerical methods for solving the NSEs have been investigated by many authors [7,8,9,10,11,12,23]. The study of time-fractional Navier-Stokes equations (TFNSEs) has become a hot topic of research due to its significant role in simulating the anomalous diffusion in fractal media. There are also some analytical methods available for solving the TFNSEs. Momani and Odibat [13] applied Adomian decomposition method to obtain the analytical solution of the TFNSEs. In [14, 15], the homotopy perturbation (transform) method was used to find the analytical solution of the TFNSEs. Wang and Liu [16] solved TFNSEs by applying the transform methods. Concerning the existence of global and local mild solutions to TFNSEs, see Carvalho-Neto and Gabriela [17], Zhou and Peng [18]. Moreover, Zhou and Peng [19] investigated the existence of weak solutions and optimal control for TFNSEs, while Peng et al.[20] presented the rigorous exposition of local solutions of TFNSEs in Sobolev space. However, one can notice that there are only a few works related to the numerical solution of the TFNSEs. The details of meshless local Petrov-Galerkin method based on moving Kriging interpolation for solving the TFNSEs can be found in the literature [21]. The purpose of this paper is to present finite difference/element method to obtain the numerical solution of TFNSEs.

The rest of the paper is arranged as follows. In Section 2, we give some notations and preliminaries. Section 3 deals with a semi-discrete scheme for the TFNSEs, which is based on a mixed finite element method in space. We also discuss the stability and error estimates of this semi-discrete scheme. In Section 4, we use a finite difference approximation to discrete time direction to obtain the fully discrete scheme. The stability and error estimates for the discrete schemes are also found. In Section 5, numerical results are discussed to confirm our theoretical analysis. Conclusions are given in the final section.

## 2. Notations and preliminaries

In this section, we present some preliminary concepts of the functional spaces. Firstly, we introduce the following Hilbert spaces:

 X=H10(Ω)2, Y=L2(Ω)2, M=L20(Ω)={v∈L2(Ω);∫Ωvdx=0},

where the space is associated with the usual inner product and the norm . The space is associated with the following inner product and equivalent norm:

 ((u,v))=(∇u,∇v), ∥u∥X=∥∇u∥0=∥u∥1.

Denote by and the closed subsets of and respectively, which are given by

 V={v∈X;div v=0}, H={v∈Y;div v=0,v⋅n|∂Ω=0}.

We denote the Stokes operator by , in which is the -orthogonal projection of onto . The domain of is and let with the norm . Observe that , and .

Next, we define the Riemann-Liouville fractional integral operator of order () as (see [22])

 Iβtg(t)=1Γ(β)∫t0(t−s)β−1g(s)ds,t>0, (2.1)

with .

The Caputo-type derivative of order , in (1.1) is defined by

 CDαtg(t)=ddt{I1−αt[g(t)−g(0)]}=ddt{1Γ(1−α)∫t0(t−s)−α[g(t)−g(0)]ds}. (2.2)

Further, the operator is defined as

 CD−αtg(t)=Iαtg(t)=1Γ(α)∫t0(t−s)α−1g(s)ds,t>0. (2.3)

where stands for the gamma function .

Next we introduce the following continuous bilinear forms and on and respectively as follows:

 a(u,v)=ν(∇u,∇v), u,v∈X, d(v,q)=(q,divv), v∈X,q∈M,

and the trilinear form on is given by

 b(u,v,w)=((u⋅∇)v+12(divu)v,w)=12((u⋅∇)v,w)−12((u⋅∇)w,v), u,v,w∈X.

It is well-known that the trilinear form has the following properties:

 b(u,v,w)=−b(u,w,v), b(u,v,v)=0, u,v,w∈X, |b(u,v,w)|≤μ0∥u∥1∥v∥1∥w∥1, u,v,w∈X.

In terms of the above notations, the weak formulation of problem (1.1) is as follows: find for all such that for all :

 {CDαt(u,v)+a(u,v)+b(u,u,v)−d(v,p)+d(u,q)=(f,v),u(0)=u0. (2.4)

In [19], Zhou and Peng discussed the existence and uniqueness of weak solutions for the problem (2.4). The objective of the present work is to obtain the numerical solution of the problem at hand.

## 3. Finite element method for space discretization

Let be a mesh of with a mesh size function , which is the diameter of element containing . Assuming be the largest mesh size of , we introduce the mixed finite element subspace of and define the subspace of as

 Vh={vh∈Xh;d(vh,qh)=0, ∀qh∈Mh}.

Let denote the -orthogonal projection defined by

 (Phv,vh)=(v,vh), ∀v∈Y, vh∈Vh}.

With the above notations, we need some further basic assumptions on the mixed finite element spaces (Refs.[8,10,12]).

(A1) Approximation. For each , there exist approximations such that

 ∥v−πhv∥1≤Ch∥Av∥, ∥q−ρhq∥≤Ch∥q∥1. (3.1)

(A2) Inverse estimate. For any , the following relations hold:

 ∥∇v∥≤Ch−1∥v∥, ∥q∥≤Ch−1∥q∥−1. (3.2)

(A3) Stability property. For any , the well-known inf-sup condition holds:

 supv∈Xhd(v,q)∥v∥1≥λ∥q∥, (3.3)

where is a constant.

Further, the following classical properties hold:

 ∥v−Phv∥+h∥∇(v−Phv)∥≤Ch2∥Av∥,v∈D(A); (3.4) ∥v−Phv∥≤Ch∥∇(v−Phv)∥,v∈X. (3.5)

The standard finite element Galerkin approximation for (2.4) holds as follows: Find for all such that for all , we have

 {CDαt(uh,vh)+a(uh,vh)+b(uh,uh,vh)−d(vh,ph)+d(uh,qh)=(f,vh),uh(0)=u0h=Phu0. (3.6)

With the above semi-discrete approximation, a discrete analogue of the Stokes operator is defined as via the condition for all . The trilinear form satisfies the following properties:

 b(uh,vh,wh)=−b(uh,wh,vh), b(uh,vh,vh)=0, uh,vh,wh∈Xh; (3.7) |b(uh,vh,wh)|≤μ0∥uh∥1∥vh∥1∥wh∥1, uh,vh,wh∈Xh. (3.8)

Theorem 3.1. For any and , let be the solution of equation (3.6). Then there exists a positive constant such that

 ∥uh∥2+ν1∫t0∥uh∥21ds≤∥u0h∥2+(1−α1)T1+β2ν(1+β)Γ(α)+α1T2νΓ(α)maxt∈[0,T]∥f∥2α1−1,

where is constant.

Proof. Taking , in (3.6) and using the Young’s inequality, we get

 CDαt∥uh∥2+ν∥uh∥21=(f,uh)≤∥f∥−1∥uh∥1≤12ν∥f∥2−1+ν2∥uh∥21, (3.9)

where denotes the dual operator in .

Applying the integral operator (2.3) to both sides of (3.9) and using the Young’s inequality, we obtain

 ∥uh∥2+ν2Γ(α)∫t0(t−s)α−1∥uh∥21ds≤∥u0h∥2+12νΓ(α)∫t0(t−s)α−1∥f∥2−1ds ≤∥u0h∥2+1−α12νΓ(α)∫t0(t−s)α−11−α1ds+α12νΓ(α)∫t0∥f∥2α1−1ds ≤∥u0h∥2+(1−α1)T1+β2ν(1+β)Γ(α)+α1T2νΓ(α)maxt∈[0,T]∥f∥2α1−1,

where with .

In view of the inequality

 ν2Γ(α)∫t0(t−s)α−1∥uh∥21ds≥νTα−12Γ(α)∫t0∥uh∥21ds,

it follows that

 ∥uh∥2+νTα−12Γ(α)∫t0∥uh∥21ds≤∥u0h∥2+(1−α1)T1+β2ν(1+β)Γ(α)+α1T2νΓ(α)maxt∈[0,T]∥f∥2α1−1.

This completes the proof.

Theorem 3.2. For any , , let and be the solutions of equations (2.4) and (3.6) respectively. Then there exists a positive constant such that

 ∥u−uh∥≤Ch2, ∥p−ph∥≤Ch. (3.10)

Proof. Setting , we deduce from (2.4) and (3.6) that

 {CDαt(ξ,v)+a(ξ,v)+b(ξ,uh,v)+b(uh,ξ,v)+b(ξ,ξ,v)−d(v,η)+d(ξ,q)=0,ξ0=u0−Phu0. (3.11)

Taking and in (3.11), we get

 CDαt∥ξ∥2+ν∥ξ∥21+b(ξ,uh,ξ)=0.

Using the properties of together with Young’s inequality, we obtain

 CDαt∥ξ∥2+ν∥ξ∥21 =b(ξ,ξ,uh) ≤C0∥ξ∥∥ξ∥1∥uh∥1 ≤ν∥ξ∥21+C1∥ξ∥2∥uh∥21. (3.12)

Applying (2.3) to both sides of (3.12), we have

 ∥ξ∥2≤∥ξ0∥2+C1Γ(α)∫t0(t−s)α−1∥ξ∥2∥uh∥21ds. (3.13)

By means of the generalized integral version of Gronwall’s lemma [24], we get

 ∥ξ∥2 ≤∥ξ0∥2exp(C1Γ(α)∫t0(t−s)α−1∥uh∥21ds) ≤∥u0−Phu0∥2exp[C2(∥u0h∥2+(1−α1)T1+β2ν(1+β)Γ(α)+α1T2νΓ(α)maxt∈[0,T]∥f∥2α1−1)] ≤Ch4. (3.14)

Furthermore, setting and in (3.11), and using inf-sup condition (3.3), combining (3.1)-(3.5),(3.8),(3.14) and using the integral operator (2.3) in (3.15), we conclude that

 ∥η∥ ≤supVh|d(ξ,η)|λ∥ξ∥1 =supVh|CDαt∥ξ∥2+ν∥ξ∥21+b(ξ,uh,ξ)|λ∥ξ∥1 ≤Ch. (3.15)

This completes the proof.

## 4. Finite difference method for time discretization

The discretization of time-fractional derivative can be found in [25-32] and references therein. Here, we will introduce a uniform grid by discretizing the temporal domain given by the points: for , with the time-step size . Hence, the Riemann-Liouville fractional integral operator of order can be discretized as follows:

 Iαtg(tn) =1Γ(α)n∑k=1∫tktk−1(tn−s)α−1g(s)ds =1Γ(α)n∑k=1(∫tktk−1(tn−s)α−1g(tk)ds)+γnα =ταΓ(α+1)n∑k=1g(tk)[(n−k+1)α−(n−k)α]+γnα =ταΓ(α+1)n−1∑k=0wαkg(tn−k)+γnα, (4.1)

where and the truncation error is given by

 γnα =1Γ(α)n∑k=1∫tktk−1(tn−s)α−1[g(s)−g(tk)]ds =1Γ(α)n∑k=1∫tktk−1(tn−s)α−1g′(ζ)(s−tk)ds, s<ζ

Therefore, we have

 |γnα| ≤τΓ(α)max0≤t≤tk|g′(t)|n∑k=1∫tktk−1(tn−s)α−1ds ≤Nατα+1Γ(α+1)max0≤t≤tk|g′(t)|.

Lemma 4.1. (see [25]) If , then

 Iαtg(tn)=ταΓ(α+1)n−1∑k=0wαkg(tn−k)+γnα, (4.2)

where .

Lemma 4.2. (see [25]) For and , let the coefficient be given by (4.1). Then

 (i) wα0=1,wαk>0,k=0,1,2,⋯; (ii) wαk>wαk+1,k=0,1,2,⋯; (iii)n−1∑k=0wαk=nα≤Nα.

Applying the integral operator (2.3) to both sides of (3.6), we obtain

 (uh,vh)+1Γ(α)∫t0(t−s)α−1[a(uh,vh)+b(uh,uh,vh)−d(vh,ph)+d(uh,qh)]ds =(u0h,vh)+1Γ(α)∫t0(t−s)α−1(f,vh)ds. (4.3)

Let and be the numerical solutions of and at respectively. By (4.2) and (4.3), our full discrete scheme of equation (2.4) can be defined by seeking such that for all :

 (unh,vh)+β0n−1∑k=0wαk[a(un−kh,vh)+b(un−kh,un−kh,vh)−d(vh,pn−kh)+d(un−kh,qh)] =(u0h,vh)+β0n−1∑k=0wαk(fn−k,vh), (4.4)

where .

Theorem 4.1. For any , the full discrete scheme (4.4) is unconditionally stable, and that

 ∥unh∥2+β1∥unh∥21≤C†(∥u0h∥2+n∑k=0∥fk∥2−1).

Proof. Setting in (4.4), we get

 (u1h,vh)+β0[a(u1h,vh)+b(u1h,u1h,vh)−d(vh,p1h)+d(u1h,qh)]=(u0h,vh)+β0(f1,vh). (4.5)

Taking and in (4.5), we have

 (u1h,u1h)+β0a(u1h,u1h)=(u0h,u1h)+β0(f1,u1h).

Making use of Young’s inequality, we obtain

 ∥u1h∥2+β0ν∥u1h∥21≤12[∥u0h∥2+∥u1h∥2]+[β02ν∥f1∥2−1+β0ν2∥u1h∥21],

that is,

 ∥u1h∥2+β1∥u1h∥21≤12∥u0h∥2+β02ν∥f1∥2−1.

Assuming and , the following inequality holds

 ∥ujh∥2+β1∥ujh∥21≤C(∥u0h∥2+j∑k=1∥fk∥2−1), j=2,3,…,n−1.

Setting and in (4.4), we get

 (unh,un−kh)+β0n−1∑k=0wαka(un−kh,un−kh)=(u0h,un−kh)+β0n−1∑k=0wαk(fn−k,un−kh).

By the elementary identity and the Young’s inequality, we have

 12[∥unh∥2+∥un−kh∥2]+β0νn−1∑k=0wαk∥un−kh∥21 =(u0h,un−kh)+12∥unh−un−kh∥2+β0n−1∑k=0wαk(fn−k,un−kh) ≤12[∥u0h∥2+∥un−kh∥2]+12∥unh−un−kh∥2+n−1∑k=0wαk(β02ν∥fn−k∥2−1+β0ν2∥un−kh∥21).

Together with Lemma 4.2 (ii) (that is, ) and , we obtain

 ∥unh∥2+β1∥unh∥21 ≤∥u0h∥2+∥unh−un−kh∥2+β02νn−1∑k=0wαk∥fn−k∥2−1 ≤C†(∥u0h∥+n∑k=0∥fk∥2−1). (4.6)

The proof is completed.

Lemma 4.3. Let be the viscosity coefficient and that

 ∥u0h∥2+n∑k=0∥fk∥2−1≤β1νC†μ0. (4.7)

Then

 a(vh,vh)+b(vh,unh,vh)≥0, (4.8)

where is defined by (3.8) and is constant.

Proof. Making use of (4.6) and (4.7), we get

 ∥unh∥21≤νμ0,

that is,

 ν−μ0∥unh∥21≥0.

By the property of , we get

 a(vh,vh)+b(vh,unh,vh)≥(ν−μ0∥unh∥21)∥vh∥21≥0.

The proof of the lemma is completed.

Theorem 4.2. For , let and be the solutions of equations (3.6) and (4.4) respectively. There exists a constant such that

 ∥uh(tn)−unh∥≤Cτα+1, ∥ph(tn)−pnh∥≤Cτα+1. (4.8)

Proof. Let and . Using (4.2)-(4.4) and noting , we deduce

 (ξn,vh)+β0n−1∑k=0wαk[a(ξn−k,vh)+b(ξn−k,un−kh,vh)−d(vh,ηn−k)+d(ξn−k,qh)] =(γnα,vh). (4.9)

For , taking and in (4.9), we have

 (ξ1,ξ1)+β0[a(ξ1,ξ1)+b(ξ1,u1h,ξ1)]=(γ1α,ξ1).

By Cauchy-Schwarz inequality and Lemma 4.3, we get

 ∥ξ1∥≤∥γ1α∥≤Cτα+1. (4.10)

Let us assume that for . In order to show that the first inequality in (4.8) holds for