Finite difference/element method for time-fractional Navier-Stokes equations
Guang-an Zou Yong Zhou Bashir Ahmad Ahmed Alsaedi







G. Zhou, Y. Zhou, B. Ahmad, a. Alsaedi: Finite difference/element method ... ] Finite difference/element method for time-fractional Navier-Stokes equations
Appl. Comput. Math., V.xx, N.xx,  20xx]Guang-an Zou^1,  Yong Zhou^2,3, Bashir Ahmad^3, Ahmed Alsaedi^3
^1School of Mathematics and Statistics, Henan University, Kaifeng 475004, P. R. China

^2Faculty of Mathematics and Computational Science, Xiangtan University
Xiangtan,  Hunan 411105, P.R. China

^3Nonlinear Analysis and Applied Mathematics (NAAM) Research Group, Faculty of Science

King Abdulaziz University, Jeddah 21589, Saudi Arabia

   e-mail: zouguangan@henu.edu.cn, yzhou@xtu.edu.cn, bashirahmad_qau@yahoo.com, aalsaedi@hotmail.com

    Manuscript received xx


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 0<α<1. 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

[
    [
    December 30, 2023
=====================





§ INTRODUCTION
In this paper, we study the following Navier-Stokes equations with time-fractional derivative in a bounded subset of Ω⊂𝐑^2 with a smooth boundary ∂Ω:

    ^CD_t^αu+(u·∇)u-νΔ u+∇ p=f,   div u=0, ∀ (x,t)∈Ω×(0,T],
     u(x,t)|_∂Ω=0, t∈ (0,T],
     u(x,0)=u_0, x∈Ω,
    1.1

where ^CD_t^α represents the Caputo-type fractional derivative of order α∈(0,1), u=(u_1(x,t), u_2(x,t)) denotes the velocity field at a point x∈Ω and t∈[0,T], ν>0 is viscosity coefficient, p=p(x,t) represents the pressure field, f=f(x,t) is the external force and u_0=u_0(x) is the initial velocity.

Notice that the problem (1.1) reduces to the classical Navier-Stokes equations (NSEs) for α=1. 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.


§ NOTATIONS AND PRELIMINARIES


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

    X=H_0^1(Ω)^2, Y=L^2(Ω)^2, M=L_0^2(Ω)={v∈ L^2(Ω);∫_Ωv dx=0},

where the space L^2(Ω) is associated with the usual inner product (·,·) and the norm ·. The space X is associated with the following inner product and equivalent norm:

    ((u,v))=(∇ u,∇ v), u_X=∇ u_0= u_1.


Denote by V and H the closed subsets of X and Y 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 A=-P, in which P is the L^2-orthogonal projection of Y onto H. The domain of A is 𝒟(A)=H^2(Ω)^2∩ V and let H^s=𝒟(A^s/2) with the norm v_s=A^s/2 v. Observe that H^2=𝒟(A), H^1=V and H^0=H.

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

    I_t^βg(t)=1/Γ(β)∫_0^t(t-s)^β-1g(s)ds,t>0, 2.1

with I_t^0g(t)=g(t).

The Caputo-type derivative of order α∈ (0,1], ^CD_t^α in (1.1) is defined by

    ^CD_t^αg(t)=d/dt{I_t^1-α[g(t)-g(0)]}=d/dt{1/Γ(1-α)∫_0^t(t-s)^-α[g(t)-g(0)]ds}.2.2

Further, the operator ^CD_t^-α is defined as

    ^CD_t^-αg(t)=I_t^αg(t)=1/Γ(α)∫_0^t(t-s)^α-1g(s)ds,t>0.2.3

where Γ(·) stands for the gamma function Γ(x)=∫_0^∞t^x-1e^-tdt.

Next we introduce the following continuous bilinear forms a(·,·) and d(·,·) on X× X and X× M respectively as follows:

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

and the trilinear form b(u,v,w) on X× X× X is given by

    b(u,v,w)=((u·∇)v+1/2(divu)v,w)=1/2((u·∇)v,w)-1/2((u·∇)w,v), u,v,w∈ X.


It is well-known that the trilinear form b(u,v,w) 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)|≤μ_0u_1v_1w_1, u,v,w∈ X.


In terms of the above notations, the weak formulation of problem (1.1) is as follows: find (u,p)∈(X,M) for all t∈ [0,T] such that for all (v,q)∈(X,M):

    ^CD_t^α(u,v)+a(u,v)+b(u,u,v)-d(v,p)+d(u,q)=(f,v),
    
    u(0)=u_0.
    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.


§ FINITE ELEMENT METHOD FOR SPACE DISCRETIZATION


Let T^h(Ω)={K} be a mesh of Ω with a mesh size function h(x), which is the diameter h_K of element K containing x.
Assuming h=h_Ω=max_x∈Ωh(x) be the largest mesh size of T^h(Ω), we introduce the mixed finite element subspace (X_h,M_h) of (X,M) and define the subspace V_h of X_h as

    V_h={v_h∈ X_h;d(v_h,q_h)=0, ∀ q_h∈ M_h}.


Let P_h:Y→ V_h denote the L^2-orthogonal projection defined by

    (P_hv,v_h)=(v,v_h), ∀ v∈ Y, v_h∈ V_h}.


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

(A1) Approximation. For each (v,q)∈ (𝒟(A),M∩ H^1(Ω)), there exist approximations (π_hv,ρ_hq)∈ (X_h,M_h) such that

    v-π_hv_1≤ ChA v, q-ρ_hq≤ Chq_1.3.1


(A2) Inverse estimate. For any (v,q)∈ (X_h,M_h), the following relations hold:

    ∇ v≤ Ch^-1 v, q≤ Ch^-1q_-1.3.2


(A3) Stability property. For any (v,q)∈ (X_h,M_h), the well-known inf-sup condition holds:

    sup_v∈ X_hd(v,q)/v_1≥λq,3.3

where λ>0 is a constant.

Further, the following classical properties hold:

    v-P_hv+h∇ (v-P_hv)≤ Ch^2A v,v∈𝒟(A);3.4
       v-P_hv≤ Ch∇ (v-P_hv),v∈ X.3.5


The standard finite element Galerkin approximation for (2.4) holds as follows: Find (u_h,p_h)∈(X_h,M_h) for all t∈ [0,T] such that for all (v_h,q_h)∈(X_h,M_h), we have

    ^CD_t^α(u_h,v_h)+a(u_h,v_h)+b(u_h,u_h,v_h)-d(v_h,p_h)+d(u_h,q_h)=(f,v_h),
    
    u_h(0)=u_0h=P_hu_0.
    3.6


With the above semi-discrete approximation, a discrete analogue of the Stokes operator A is defined as A_h=-P_h_h via the condition (-_hu_h,v_h)=((u_h,v_h)) for all u_h,v_h∈ X_h. The trilinear form b(u_h,v_h,w_h) satisfies the following properties:

    b(u_h,v_h,w_h)=-b(u_h,w_h,v_h), b(u_h,v_h,v_h)=0, u_h,v_h,w_h∈ X_h;3.7
       |b(u_h,v_h,w_h)|≤μ_0u_h_1v_h_1w_h_1, u_h,v_h,w_h∈ X_h.3.8


Theorem 3.1. For any t∈ [0,T] and 0<α < 1, let u_h be the solution of equation (3.6). Then there exists a positive constant C such that

    u_h^2+ν_1∫_0^tu_h_1^2ds≤u_0h^2+(1-α_1)T^1+β/2ν(1+β)Γ(α)+α_1T/2νΓ(α)max_t∈ [0,T]f_-1^2/α_1,

where ν_1=ν T^α-1/2Γ(α)>0 is constant.

Proof. Taking v_h=u_h, q_h=p_h in (3.6) and  using the Young's inequality, we get

    ^CD_t^αu_h^2+νu_h_1^2=(f,u_h) ≤f_-1u_h_1≤1/2νf_-1^2+ν/2u_h_1^2,3.9

where f_-1=A^-1/2f denotes the dual operator in 𝒟(A^-1/2).

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

    u_h^2+ν/2Γ(α)∫_0^t(t-s)^α-1u_h_1^2ds≤u_0h^2+1/2νΓ(α)∫_0^t(t-s)^α-1f_-1^2ds
       ≤u_0h^2+1-α_1/2νΓ(α)∫_0^t(t-s)^α-1/1-α_1ds+α_1/2νΓ(α)∫_0^tf_-1^2/α_1ds
       ≤u_0h^2+(1-α_1)T^1+β/2ν(1+β)Γ(α)+α_1T/2νΓ(α)max_t∈ [0,T]f_-1^2/α_1,

where β=α-1/1-α_1 with 0<α_1 < 1.

In view of the inequality

    ν/2Γ(α)∫_0^t(t-s)^α-1u_h_1^2ds≥ν T^α-1/2Γ(α)∫_0^tu_h_1^2ds,

it follows that

    u_h^2+ν T^α-1/2Γ(α)∫_0^tu_h_1^2ds≤u_0h^2+(1-α_1)T^1+β/2ν(1+β)Γ(α)+α_1T/2νΓ(α)max_t∈ [0,T]f_-1^2/α_1.


This completes the proof.

Theorem 3.2. For any t∈ [0,T], 0<α < 1, let (u,p) and (u_h,p_h) be the solutions of equations (2.4) and (3.6) respectively. Then there exists a positive constant C such that

    u-u_h≤ Ch^2, p-p_h≤ Ch. 3.10



Proof. Setting (ξ,η)=(u-u_h,p-p_h), we deduce from (2.4) and (3.6) that

    ^CD_t^α(ξ,v)+a(ξ,v)+b(ξ,u_h,v)+b(u_h,ξ,v)+b(ξ,ξ,v)-d(v,η)+d(ξ,q)=0,
    ξ_0=u_0-P_hu_0.
    3.11


Taking v=ξ and q=η in (3.11), we get

    ^CD_t^αξ^2+νξ_1^2+b(ξ,u_h,ξ)=0.


Using the properties of b(u_h,v_h,w_h) together with Young's inequality, we obtain

    ^CD_t^αξ^2+νξ_1^2   =b(ξ,ξ,u_h)
       ≤ C_0ξξ_1u_h_1
       ≤νξ_1^2+C_1ξ^2u_h_1^2.3.12


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

    ξ^2≤ξ_0^2+C_1/Γ(α)∫_0^t(t-s)^α-1ξ^2u_h_1^2ds.3.13


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

    ξ^2   ≤ξ_0^2exp(C_1/Γ(α)∫_0^t(t-s)^α-1u_h_1^2ds)
       ≤u_0-P_hu_0^2exp[C_2(u_0h^2+(1-α_1)T^1+β/2ν(1+β)Γ(α)+α_1T/2νΓ(α)max_t∈ [0,T]f_-1^2/α_1)]
       ≤ Ch^4.3.14



Furthermore, setting v=ξ and q=0 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

    η   ≤sup_V_h|d(ξ,η)|/λξ_1
       =sup_V_h|^CD_t^αξ^2+νξ_1^2+b(ξ,u_h,ξ)|/λξ_1
       ≤ Ch. 3.15


This completes the proof.



§ 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 [0,T] given by the points: t_n=nτ for n=0,1,…,N, with the time-step size τ=T/N. Hence, the Riemann-Liouville fractional integral operator of order α can be discretized as follows:

    I_t^αg(t_n)   =1/Γ(α)∑_k=1^n∫_t_k-1^t_k(t_n-s)^α-1g(s)ds
       =1/Γ(α)∑_k=1^n(∫_t_k-1^t_k(t_n-s)^α-1g(t_k)ds)+γ_α^n
       =τ^α/Γ(α+1)∑_k=1^ng(t_k)[(n-k+1)^α-(n-k)^α]+γ_α^n
       =τ^α/Γ(α+1)∑_k=0^n-1w_k^αg(t_n-k)+γ_α^n,4.1

where w_k^α=(k+1)^α-k^α and the truncation error γ_α^n is given by

    γ_α^n   =1/Γ(α)∑_k=1^n∫_t_k-1^t_k(t_n-s)^α-1[g(s)-g(t_k)]ds
       =1/Γ(α)∑_k=1^n∫_t_k-1^t_k(t_n-s)^α-1g^'(ζ)(s-t_k)ds, s<ζ<t_k.


Therefore, we have

    |γ_α^n|   ≤τ/Γ(α)max_0≤ t≤ t_k|g^'(t)|∑_k=1^n∫_t_k-1^t_k(t_n-s)^α-1ds
       ≤N^ατ^α+1/Γ(α+1)max_0≤ t≤ t_k|g^'(t)|.


Lemma 4.1. (see <cit.>) If g(t)∈ C^1[0,T], then

    I_t^αg(t_n)=τ^α/Γ(α+1)∑_k=0^n-1w_k^αg(t_n-k)+γ_α^n,  4.2

where |γ_α^n|≤ Cτ^α+1, n=0,1,…,N.

Lemma 4.2. (see <cit.>) For 0< t_n≤ t_N=T and α>0, let the coefficient w_k^α 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)∑_k=0^n-1w_k^α=n^α≤ N^α.



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

    (u_h,v_h)+1/Γ(α)∫_0^t(t-s)^α-1[a(u_h,v_h)+b(u_h,u_h,v_h)-d(v_h,p_h)+d(u_h,q_h)]ds
       =(u_0h,v_h)+1/Γ(α)∫_0^t(t-s)^α-1(f,v_h)ds.  4.3


Let u_h^n and p_h^n be the numerical solutions of u_h(t) and p_h(t) at t=t_n respectively. By (4.2) and (4.3), our full discrete scheme of equation (2.4) can be defined by seeking
(u_h^n,p_h^n)∈(X_h,M_h) such that for all (v_h,q_h)∈(X_h,M_h):

    (u_h^n,v_h)+β_0∑_k=0^n-1w_k^α[a(u_h^n-k,v_h)+b(u_h^n-k,u_h^n-k,v_h)-d(v_h,p_h^n-k)+d(u_h^n-k,q_h)]
       =(u_h^0,v_h)+β_0∑_k=0^n-1w_k^α(f^n-k,v_h), 4.4

where β_0=τ^α/Γ(α+1).


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

    u_h^n^2+β_1u_h^n_1^2≤ C^†(u_h^0^2+∑_k=0^nf^k_-1^2).


Proof. Setting n=1 in (4.4), we get

    (u_h^1,v_h)+β_0[a(u_h^1,v_h)+b(u_h^1,u_h^1,v_h)-d(v_h,p_h^1)+d(u_h^1,q_h)]=(u_h^0,v_h)+β_0(f^1,v_h). 4.5


Taking v_h=u_h^1 and q_h=p_h^1 in (4.5), we have

    (u_h^1,u_h^1)+β_0a(u_h^1,u_h^1)=(u_h^0,u_h^1)+β_0(f^1,u_h^1).


Making use of Young's inequality, we obtain

    u_h^1^2+β_0νu_h^1_1^2≤1/2[u_h^0^2+u_h^1^2]+[β_0/2νf^1_-1^2+β_0ν/2u_h^1_1^2],

that is,

    u_h^1^2+β_1u_h^1_1^2≤1/2u_h^0^2+β_0/2νf^1_-1^2.


Assuming v_h=u_h^j∈ X_h and q_h=p_h^j∈ M_h, the following inequality holds

    u_h^j^2+β_1u_h^j_1^2≤ C(u_h^0^2+∑_k=1^jf^k_-1^2), j=2,3,…,n-1.



Setting v_h=u_h^n-k and q_h=p_h^n-k in (4.4), we get

    (u_h^n,u_h^n-k)+β_0∑_k=0^n-1w_k^αa(u_h^n-k,u_h^n-k)=(u_h^0,u_h^n-k)+β_0∑_k=0^n-1w_k^α(f^n-k,u_h^n-k).


By the elementary identity ab=1/2(a^2+b^2)-1/2(a-b)^2 and the Young's inequality, we have

    1/2[u_h^n^2+u_h^n-k^2]+β_0ν∑_k=0^n-1w_k^αu_h^n-k_1^2
       =(u_h^0,u_h^n-k)+1/2u_h^n-u_h^n-k^2
    +β_0∑_k=0^n-1w_k^α(f^n-k,u_h^n-k)
       ≤1/2[u_h^0^2+u_h^n-k^2]+1/2u_h^n-u_h^n-k^2+∑_k=0^n-1w_k^α(β_0/2νf^n-k_-1^2+β_0ν/2u_h^n-k_1^2).


Together with Lemma 4.2 (ii) (that is, 0<w_k+1^α<w_k^α<1) and β_0ν/2∑_k=1^n-1w_k^αu_h^n-k_1^2≥0, we obtain

    u_h^n^2+β_1u_h^n_1^2   ≤u_h^0^2+u_h^n-u_h^n-k^2
    +β_0/2ν∑_k=0^n-1w_k^αf^n-k_-1^2
       ≤ C^†(u_h^0+∑_k=0^nf^k_-1^2). 4.6


The proof is completed.

Lemma 4.3. Let ν>0 be the viscosity coefficient and that

    u_h^0^2+∑_k=0^nf^k_-1^2≤β_1ν/C^†μ_0. 4.7

Then

    a(v_h,v_h)+b(v_h,u_h^n,v_h)≥0, 4.8

where μ_0>0 is defined by (3.8) and C^†>0 is constant.


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

    u_h^n_1^2≤ν/μ_0,

that is,

    ν-μ_0u_h^n_1^2≥0.


By the property of b(v_h,u_h^n,v_h), we get

    a(v_h,v_h)+b(v_h,u_h^n,v_h)≥ (ν-μ_0u_h^n_1^2)v_h_1^2≥ 0.


The proof of the lemma is completed.

Theorem 4.2. For 0<α < 1, let (u_h(t_n),p_h(t_n)) and (u_h^n,p_h^n) be the solutions of equations (3.6) and (4.4) respectively. There exists a constant C such that

    u_h(t_n)-u_h^n≤ C τ^α+1, p_h(t_n)-p_h^n≤ C τ^α+1. 4.8


Proof. Let ξ^n=u_h(t_n)-u_h^n and η^n=p_h(t_n)-p_h^n. Using (4.2)-(4.4) and noting ξ^0=0, we deduce

    (ξ^n,v_h)+β_0∑_k=0^n-1w_k^α[a(ξ^n-k,v_h)+b(ξ^n-k,u_h^n-k,v_h)-d(v_h,η^n-k)+d(ξ^n-k,q_h)]
       =(γ_α^n,v_h). 4.9


For n=1, taking v_h=ξ^1 and q_h=η^1 in (4.9), we have

    (ξ^1,ξ^1)+β_0[a(ξ^1,ξ^1)+b(ξ^1,u_h^1,ξ^1)]=(γ_α^1,ξ^1).


By Cauchy-Schwarz inequality and Lemma 4.3, we get

    ξ^1≤γ_α^1≤ C τ^α+1.  4.10


Let us assume that ξ^m≤ C τ^α+1 for m=2,3,…,n-1. In order to show that the first inequality in (4.8) holds for m=n, we set v_h=ξ^n-k and q_h=η^n-k in (4.9). Then

    (ξ^n,ξ^n-k)+β_0∑_k=0^n-1w_k^α[a(ξ^n-k,ξ^n-k)+b(ξ^n-k,u_h^n-k,ξ^n-k)]=(γ_α^n,ξ^n-k).


In view of the elementary identity ab=1/2(a^2+b^2)-1/2(a-b)^2, Young's inequality and Lemma 4.3, we get

    1/2[ξ^n^2+ξ^n-k^2]≤1/2ξ^n-ξ^n-k^2+1/2[ξ^n-k^2+γ_α^n^2], 4.11

which implies that

    ξ^n^2≤ C τ^2α+2. 4.12


By inverse estimate (3.2) together  with (4.12), we have

    ξ^n_1≤ Ch^-1ξ^n≤ C τ^α+1. 4.13



On the other hand, setting v_h=ξ^1 and q_h=0 for n=1 in (4.9) and  making use of Cauchy-Schwarz inequality together with (3.1)-(3.4), (3.8) and (4.13), we get

    η^1   ≤sup_V_h|d(ξ^1,η^1)|/λξ^1_1=sup_V_h|ξ^1^2+β_0[νξ^1_1^2+b(ξ^1,u_h^1,ξ^1)]-(γ_α^1,ξ^1)|/λξ^1_1
       ≤ C τ^α+1. 4.14


Using the assumption  η^m≤ C τ^α+1 for m=2,3,…,n-1 and taking v_h=ξ^n-k and q_h=0 in (4.9), by (4.12) and (4.13), similar to the derivation of (4.14) for m=n, we obtain

    η^n   ≤sup_V_h|d(ξ^n,η^n)|/λξ^n_1≤sup_V_h|∑_k=0^n-1w_k^αd(ξ^n-k,η^n-k)|/λξ^n_1
       =sup_V_h|(ξ^n,ξ^n-k)+β_0∑_k=0^n-1w_k^α[νξ^n-k_1^2+b(ξ^n-k,u_h^n-k,ξ^n-k)]-(γ_α^n,ξ^n-k)|/λξ^n_1
       ≤ C τ^α+1.



This completes the proof.

Next we give the error estimate for fully discrete scheme.

Theorem 4.3. For 0<α < 1, let (u(t_n),p(t_n)) and (u_h^n,p_h^n) be the solutions of equations (2.4) and (4.4) respectively. Then there exists a positive constant C such that

    u(t_n)-u_h^n≤ C(h^2+τ^α+1), p(t_n)-p_h^n≤ C(h+τ^α+1). 4.15


Proof. It is easy to show that (4.15) follows from Theorem 3.2 and Theorem 4.2 via triangle inequality.




§ NUMERICAL EXAMPLE


In this section, we demonstrate the effectiveness of our numerical methods with the aid of examples. We use mixed finite element method for the discretization of spatial direction and finite difference approximation for time discretization. The convergence rates of numerical solutions with respect to space step h and time step τ are discussed. We consider the regular (uniform) domain Ω=(0,1)×(0,1) and the time interval is chosen to be [0,1] with  the viscosity coefficient ν=1.5.
For an appropriate body force f, the analytical solution (u,p)=((u_1,u_2),p) of the unstable flow problem with homogeneous boundary conditions becomes

    u_1=2x^2(x-1)^2y(y-1)(2y-1)e^-t, u_2=-2y^2(y-1)^2x(x-1)(2x-1)e^-t,
       p=(x^2-y^2)e^-t,

which automatically satisfy the initial and boundary conditions.

The errors e^n are computed in L^2-discrete norm. The results of numerical experiments are compared with analytical solution by the rates of the convergence, which are approximately by

    Rate = |ln(e_f^n/e_c^n)/ln(N_f/N_c)|,

where e_f^n and e_c^n denote the error on finer grid and coarser grid, N_f and  N_c represent the numbers of meshes on finer grid and coarser grid, respectively.

The spatial convergence rates for the components of velocity (u_1,u_2) and pressure p with fixed time step τ=1/8 with different values of α are shown in Fig.1. The convergence rates of velocity (u_1,u_2) are in accordance with spatial convergence order 𝒪(h^2) and the pressure p are closer to order 𝒪(h). Fig.2 give the temporal convergence rates for the components of velocity and pressure with fixed spatial step h=1/15 with different values of α. We can see that the rates of convergence are closer to the theoretical convergence order 𝒪(τ^α+1).


Fig.3 depicts the numerical solutions of the components of velocity (u_1,u_2) and pressure p, with h=1/15 and τ=0.1, when α=0.4 and α=0.8, respectively. It is not difficult to find that a pair of warm- and cold-core eddies emerge in the velocity field.










§ CONCLUSION


In this study, the finite difference/element method is presented to solve the TFNSEs and the convergence error estimates for the discrete schemes in L^2-norm are obtained. We present the numerical experiment to illustrate the accuracy of schemes, and the result fully verify the convergence theory. The numerical examples also confirm the thesis [18,19,20] that in
procedure of citing and novelty of the obtained results. Furthermore, the presented methods and analytical techniques in this work can also be extended to other nonlinear time-fractional partial differential equations.








§ ACKNOWLEDGMENT

Guang-an Zou is supported by National Nature Science Foundation of China (Grant No. 11626085), Yong Zhou is supported by National Nature Science Foundation of China (Grant No. 11671339).






§ REFERENCES


[1]  Lemarié-Rieusset, P.G. Recent developments in the Navier-Stokes problem,  Chapman  Hall/CRC Research Notes in Mathematics, 431. Chapman  Hall/CRC, Boca Raton, FL, 2002, 395 p.

[2]  Chemin, J.Y.,   Gallagher, I., Paicu, M.  Global regularity for some classes of large solutions to the Navier-Stokes equations, Ann. of Math. (2), V.173, N.2, 2011, pp.983-1012.


[3] Miura, H. Remark on uniqueness of mild solutions to the Navier-Stokes equations, J. Funct. Anal., V.218, N.1, 2005, pp.110-129.

[4] Germain, P. Multipliers, paramultipliers, and weak-strong uniqueness for the Navier-Stokes equations, J. Differential Equations, V.226, N.2, 2006, pp.373-428.

[5] Nahmod, A.R., Pavlovic N., Staffilani, G. Almost sure existence of global weak solutions for supercritical Navier-Stokes equations, SIAM J. Math. Anal. V.45, N.6, 2013, pp.3431-3452.

[6] Robinson, J.C., Sadowski, W.,  Silva, R.P.  Lower bounds on blow up solutions of the three-dimensional Navier-Stokes equations in homogeneous Sobolev spaces, J. Math. Phys., V.53, N.11,  2012, 115618, 15 pp.

[7] Ingram, R. A new linearly extrapolated Crank-Nicolson time-stepping scheme for the Navier-Stokes equations, Math. Comp., V.82, N.284, 2013, pp.1953-1973.

[8] Bernardi, C.,  Raugel, G.A.  conforming finite element method for the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal., V.22, N.3, 1985, pp.455-473.

[9] He, Y., Sun, W. Stability and convergence of the Crank-Nicolson/Adams-Bashforth scheme for the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal., V.45, N.2, 2007, pp.837-869.

[10] Girault, V., Raviart, P.A. Finite element methods for Navier-Stokes equations. Theory and algorithms, Springer-Verlag, Berlin, 1986, 374 p.


[11] Kaya, S.,  Rivière, B.  A discontinuous subgrid eddy viscosity method for the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal., V.43, N.4, 2005, pp.1572-1595.

[12] Shan, L.,  Hou, Y. A fully discrete stabilized finite element method for the time-dependent Navier-Stokes equations, Appl. Math. Comput., V.215, N.1, 2009, pp.85-99.

[13] Momani, S., Odibat, Z. Analytical solution of a time-fractional Navier-Stokes equation by Adomian decomposition method, Appl. Math. Comput., V.177, N.2, 2006, pp.488-494.

[14] Ganji, Z.Z., Ganji, D.D., Ganji, Ammar D., Rostamian, M.  Analytical solution of time-fractional Navier-Stokes equation in polar coordinate by homotopy perturbation method, Numer. Methods Partial Differential Equations, V.26, N.1, 2010, pp.117-124.

[15] Kumar, D.,  Singh, J.,  Kumar, S. A fractional model of Navier-Stokes equation arising in unsteady flow of a viscous fluid, J. Ass. Arab Univ. Basic Appl. Sci., V.17, 2015, pp.14-19.

[16] Wang, K., Liu, S. Analytical study of time-fractional Navier-Stokes equation by using transform methods, Adv. Differential Equ., N.61, 2016, pp.12.

[17]  De Carvalho-Neto, P.M., Gabriela, P.  Mild solutions to the time fractional Navier-Stokes equations in R^N, J. Differential Equations, V.259, N.7, pp.2948-2980.

[18] Zhou, Y., Peng, L. On the time-fractional Navier-Stokes equations,  Comput. Math. Appl., V.73, N.6, 2017, pp.874-891.

[19] Zhou, Y., Peng, L. Weak solutions of the time-fractional Navier-Stokes equations and optimal control, Comput. Math. Appl., V.73, N.6, 2017, pp.1016-1027.

[20] Peng, L.,  Zhou, Y., Ahmad, B.,   Alsaedi, A. The Cauchy problem for fractional Navier-Stokes equations in Sobolev spaces, Chaos Solitons Fractals, V.102, 2017, pp.218-228.

[21] Thamareerat, N., Luadsong, A., Aschariyaphotha, N.  The meshless local Petrov-Galerkin method based on moving Kriging interpolation for solving the time fractional Navier-Stokes equations, SpringerPlus,  V.5, N.417, 2016, pp.19.

[22] Kilbas, A.A.,  Srivastava, H.M., Trujillo, J.J. Theory and applications of fractional differential equations, Elsevier, 2006, 523 p.

[23]  Layton, W.J., Labovschii, A.,   Manica, C.C., Neda, M., Rebholz, L. G. The stabilized, extrapolated
trapezoidal finite element method for the Navier-Stokes equations,  Comput. Methods Appl. Mech. Eng., V.198,  2009, pp.958-974.

[24] Kruse R. Strong and weak approximation of semilinear stochastic evolution equations,  Springer, Cham, 2014, 177 p.

[25] Zeng, F., Li, C., Liu, F., Turner, I. The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput., V.35, N.6, 2013, pp.A2976-A3000.

[26] Stynes, M.,  O'Riordan, E.,  Gracia, J.L.  Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal., V.55, N.2, pp.2017, 1057-1079.

[27] Stynes, M.,  Gracia, J.L. Preprocessing schemes for fractional-derivative problems to improve their convergence rates, Appl. Math. Lett., V.74, 2017, pp.187-192.

[28] Kopteva, N.,  Stynes, M. Analysis and numerical solution of a Riemann-Liouville fractional derivative two-point boundary value problem, Adv. Comput. Math.,  V.43, N.1, 2017, 77-99.

[29] Zeng, F., Li, C.,  Liu, F.,  Turner, I. Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy, SIAM J. Sci. Comput., V.37, N.1, 2015, pp.A55-A78.

[30] Zheng, M., Liu, F.,  Liu, Q.,  Burrage, K.,   Simpson, M.J. Numerical solution of the time fractional reaction-diffusion equation with a moving boundary, J. Comput. Phys.,  V.338, 2017, pp.493-510.

[31] Cui, M.R. Compact alternating direction implicit method for two-dimensional time fractional diffusion equation, J. Comput. Phys., V.231, N.6, 2012, pp.2621-2633.

[32] Jiang, Y.J.,  Ma, J.T.  High-order finite element methods for time-fractional partial differential
equations, J. Comput. Appl. Math., V.235, N.11, 2011, 3285-3290.




















































