Numerical analysis for Navier–Stokes equations with time fractional derivatives

Applied Mathematics and Computation 336 (2018) 481–489

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical analysis for Navier–Stokes equations with time fractional derivativesR Jun Zhang a, JinRong Wang b,c,∗ a

Computational mathematics research center, Guizhou University of Finance and Economics, Guiyang, Guizhou 550025, PR China Department of Mathematics, Guizhou University, Guiyang, Guizhou 550025, PR China c School of Mathematical Sciences, Qufu Normal University, Qufu 273165, PR China b

a r t i c l e

i n f o

Keywords: Navier–Stokes equations Caputo fractional derivative Finite difference Legendre-spectral method Error estimate

a b s t r a c t In this article, we study numerical approximation for a class of Navier–Stokes equations with time fractional derivatives. We propose a scheme using finite difference approach in fractional derivative and Legendre-spectral method approximations in space and prove that the scheme is unconditionally stable. In addition, the error estimate shows that the numerical solutions converge with the order O (t 2−α + t −α N 1−s ), 0 < α < 1 being the order of the fractional derivative in time. Numerical examples are illustrated to verify the theoretical results. © 2018 Elsevier Inc. All rights reserved.

1. Introduction This paper is devoted to investigate a numerical scheme for the approximation of following time-fractional Navier–Stokes equations (TFNSE):

Dtα u − ν u + u · ∇ u + ∇ p = f , in  × [0, T ],

(1.1)

∇ · u = 0, in  × [0, T ], u|t=0 = u0 ,

(1.2)

∂  × [0, T ],

(1.3)

u = 0,

on

where α ∈ (0, 1) and Dtα denotes the Caputo fractional derivative,  ⊂ R2 is a bound domain, u = (u1 , u2 ) represents the velocity, and p which represents the pressure. ν > 0 is the viscosity coefficient and f = ( f 1 , f2 ) is an external force. TFNSE have been studied extensively due to their applications in polymer physics, electrochemistry of corrosion and many other physical processes. Some recent theoretical and numerical results show that the classical diffusion equation can not describe the phenomena in heterogeneous porous media with fractal characteristics, but fractional differential equations are effective in simulating anomalous diffusion process. Fractional derivatives are found to be quite flexible in describing R This work is partially supported by National Natural Science Foundation of China (11661016) and Guizhou University of Finance and Economics (No.2017XZD01). ∗ Corresponding author at: Department of Mathematics, Guizhou University, Guiyang, Guizhou 550025, PR China. E-mail addresses: [email protected] (J. Zhang), [email protected], [email protected] (J. Wang).

https://doi.org/10.1016/j.amc.2018.04.036 0 096-30 03/© 2018 Elsevier Inc. All rights reserved.

482

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

complex dynamics in viscoelastic fluids, and the hydrodynamic model with fractional derivative can eliminate the shortage of continuous flow hypothesis [1–5]. Based on the experimental data, He [6] proposed a fractional partial differential equation for porous media percolation. Liu et al. [7–11] presented a series of numerical methods to solve the viscoelastic fluid problem with fractional derivative. The problem (1.1)–(1.3) turns into classical Navier–Stokes equations when α = 1, the related numerical methods of which have been found in [12–18] and the references therein. Recently, El-Shahed and Salem [19] obtained the analytical solution for TFNSE by using Laplace, Fourier Sine transforms, and finite Hankel transforms. Carvalhoneto and Planas [20] studied existence and uniqueness properties of local mild solutions to the TFNSE. It is remarkable that Zhou and Peng [21,22] investigated the existence and uniqueness of global and local mild solutions of TFNSE. In spite of basic theory analysis for TFNSE, some recent contributions focus on using numerical methods to approximate the solutions of TFNSE. For example, Ganji et al. [23], Kumar et al. [24] have introduced a homotopy perturbation method to simulate differential equations of Caputo fractional order arising in fluid mechanics. Kumar et al. [25] have discussed the approximated analytical solutions of TFNSE by using modified fractional Laplace decomposition method. In [26], a mixed finite element method was proposed and analyzed, under some certain assumptions on the u and f, they derived some error estimates for numerical schemes. Momani and Obibat [27] proposed a Adomian decomposition method for TFNSE in a tube. On the other hand, many numerical techniques are used to investigate fractional differential equations. Lin and Xu [28] proposed a finite difference/Legendre collocation spectral method to approximate Caputo fractional equation, and their numerical method leads to 2-α order accuracy in time and spectral accuracy in space. Li and Xu [29] constructed a timespace spectral method for fractional partial differential equations. Bhrawy et al. [30] presented a shifted Legendre collocation method to solve a time-space fractional Burgers’ equation in Caputo sense. Bhrawy et al. [31,32] introduced a Legendre spectral tau method to investigate the time-space fractional diffusion equation. At present, there is still a lack of effective numerical method to investigate TFNSE. In this article, we introduce and discuss a efficient numerical scheme to solve TFNSE. The proposed scheme is performed combining finite difference scheme for fractional derivative and spectral discretization for spatial variable. A detailed analysis of the numerical scheme is provided for both stability and error estimate. Under the appropriate assumptions, our rigorous analysis result show that the scheme is unconditionally stable, and that the convergent order is O (t 2−α + t −α N 1−s ), where t, N and s are, respectively step of time, polynomial degree, and regularity of u. At last, some numerical experiments are conducted to confirm the theoretical claims. The rest of the article is organized in following way. Section 2 introduces a scheme for fractional derivative. In Section 3, we discuss a error estimate of the full discrete scheme. In Section 4, we present some numerical experiments to illustrate the validity of the numerical method. The conclusions of this paper is given in Section 5. 2. Stability for semi-discretization fractional Navier–Stokes equations We now introduce some of the notations in this article. We use the standard notation L2 (), Hs (), and H0s () to denote the usual Sobolev spaces. We use  ·  to denote the norm in L2 (), and ( · , · ) to denote the scalar product in L2 (). The dual space of H01 () will be denote by H −1 (), and the duality between them will be denoted by ., .. We also denote the following spaces

X = {q ∈ L2 () :

 

qd = 0}.

(2.1)

And a trilinear operator B( · , · ) by

B(u, v, w ) = (u · ∇ v, w ), u, v, w ∈ H01 ().

(2.2)

Therefore, we have

B(u, v, v ) = 0, u, v, ∈ H01 ().

(2.3)

From [33] we have the following inequality:

B(u, v, w ) ≤ C0 ∇ u∇ v∇ w, ∀u, v, w ∈ H01 ().

(2.4)

In this context, we consider the discrete Eq. (1.1) and (1.2) in time. First, we introduce a scheme to discrete the Caputo derivative. For a given number T > 0 and integer M, let tn = nt, n = 0, 1, . . . , M − 1, where t = T /M. From [28,34], we know that the Caputo derivative Dtα u can be approximated by n  u(tn+1− j ) − u(tn− j ) 1 bj , (2 − α ) t α j=0

where b j = ( j + 1 )1−α − j1−α .

(2.5)

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

483

Then we obtain the following finite difference scheme for time discretization of (1.1)–(1.3):



u

n+1

−

n−1 



(b j − b j+1 )u

n− j

− bn u

− a0 ν un+1

0

(2.6)

j=0

+ a0 (un+1 · )∇ un+1 + a0 ∇ pn+1 = a0 fn+1 ,

∇ · un+1 = 0,

(2.7)

where a0 = t α (2 − α ). Lemma 2.1. The semi discrete scheme (2.6) and (2.7) is uncondition stable, for all n with 0 ≤ n ≤ M − 1,

un+1 2 + a0 ν∇ un+1 2 ≤u0 2 + Ca0

n+1 

f j 2−1 .

(2.8)

j=1

Proof. For n = 0, taking the inner product of (2.6) with 2u1 , and of (2.7) with 2a0 p1 , and making use the property of B( · , · ), we derive

( u 1 − u 0 , 2 u 1 ) + 2 a0 ν ( ∇ u 1 , ∇ u 1 ) = 2 a0 f1 , u 1 . Then a simple rearrangement yields

u1 2 − u0 2 + u1 − u0 2 + 2ν a0 ∇ u1 2 ≤ ν a0 ∇ u1 2 + Ca0 f1 2−1 , that is

u1 2 + ν a0 ∇ u1 2 ≤ u0 2 + Ca0 f1 2−1 . Suppose

uk 2 + ν a0 ∇ uk 2 ≤ u0 2 + Ca0

k 

f j 2−1 , k = 1, 2, . . . , n.

(2.9)

j=1

For k = n + 1, we have

2u



 + 2ν a0 ∇ u

n+1 2

n−1 

 =2

n+1 2

 (b j − b j+1 )u

n− j

0

+ bn u , u

n+1

+ 2a0 fn+1 , un+1 

j=0

≤

n−1 

(b j − b j+1 )(un− j 2 + un+1 2 ) + bn (u0 2 + un+1 2 ))

j=0

+ Ca0 fn+1 2−1 + ν a0 ∇ un+1 2 . Noting that n−1 

(b j − b j+1 ) + bn = 1.

j=0

By using (2.9), we get

u



 + ν a0 ∇ u

n+1 2

 ≤

n+1 2

n−1 

 u  + Ca0

(b j − b j+1 ) + bn

0 2

j=0

=u0 2 + Ca0

n 

 f(t j )

2 −1

+ Ca0 fn+1 2−1

j=1 n+1 

f j 2−1 .

j=1

The proof is completed.



Remark 2.1. If we taking the inner product of (2.6) with 2tun+1 , and of (2.7) with t pn+1 , we have following energy inequalities:

E (um+1 ) − E (u1 ) + ν t

m 

∇ un+1 2 ≤ C t

n=1

m 

fn+1 2−1 + T α u0 2 ,

n=1

where

E (um+1 ) = μ

m 

b j um+1− j 2 ,

μ = t 1−α (2 − α ).

j=0

But it is not sufficient for obtain any meaningful error estimate. Moreover, due to the discrete scheme for Caputo fractional derivative, Lemma 2.1 does not provide a k-independent stability result for ∇ pn .

484

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

3. Error estimates for full discretization In this part, we will exhibit a Legendre–Galerkin spectral method for the discretization of (2.6) and (2.7), the use of Legendre–Galerkin method is motivated by the memory feature of the fractional Navier–Stokes equation, and the Legendre– Galerkin spectral proceeds to approximate the solution by the polynomials of hight degree. First, we define SN as the polynomials space of degree less that or equal to N, let

SN0 = H01 () ∩ SN . 0 , such that We denote by π N the usual L2 -projection operator, i.e., for all v ∈ L2 (), πN v ∈ SN

( πN v − v , ψ ) = 0 ,

∀ψ ∈ SN .

0 by Similarly we define the H1 -projection operator πN1 : H01 () → SN

(∂x (πN1 v − v ), ∂x ψ ) = 0,

∀ψ ∈ SN0 .

It is known that the following estimates hold [35]:

v − πN v0 ≤cN−s vs , ∀v ∈ H0s (), s ≥ 0;

(3.1)

v − πN1 vk ≤cNk−s vs , k = 0, 1; ∀v ∈ H0s (), s ≥ 1.

(3.2)

0 , X ), such that We consider the Legendre–Galerkin spectral discretization to (2.6) and (2.7) as follows: find (unN+1 , pnN+1 ) ∈ SN 0,X) for all ( , q ) ∈ (SN

  n−1  t −α n− j n+1 0 u − (b j − b j+1 )uN − bn uN , + ν (∇ unN+1 , ∇ ), (2 − α ) N j=0

(3.3)

+B(unN+1 , unN+1 , ) + (∇ pnN+1 , ) = (fnN+1 , ),

(∇ unN+1 , q ) = 0.

(3.4)

Remark 3.1. We assume that the couple

(SN0 , X )

satisfies the discrete LBB condition.

We now state a stability result of the scheme (3.3) and (3.4). Theorem 3.1. Suppose {unN }M−1 be the solution of (3.3) and (3.4), then we have that k=1

unN+1 2 + a0 ν∇ unN+1 2 ≤u0N 2 + Ca0

n+1 

fNj 2−1 .

(3.5)

j=1

For simplification of presentation, we denoted the truncation error by

r1n+1 =

1 (1 − α )



tn+1 0

n  u(tn+1− j ) − u(tn− j ) 1 ∂ u (s ) ds − b , α ∂ s (tn+1 − s ) (2 − α ) j=0 j t α

(3.6)

r2n+1 = fn+1 − fnN+1 , r3n+1 = p(tn+1 ) − pnN+1 ,

(3.7)

r n+1 = r1n+1 + r2n+1 .

(3.8)

Following [35,36], We have the following projection estimates

r1n+1  ≤ C t 2−α , r2n+1  ≤ CN−s , r3n+1  ≤ CN−s .

(3.9)

We define the following error functions:

 enN = πN1 u(tn ) − unN ,  enN = u(tn ) − πN1 u(tn ), enN =  enN +  enN = u(tn ) − unN .

(3.10)

The next two lemmas will help us to prove the error estimate of the problem (3.3) and (3.4). Lemma 3.1. Let ν > 0, suppose that

u0N 2 + Ca0

n+1  j=1

fNj 2−1 ≤

a0 ν 3 , 4C02

(3.11)

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

485

then we have

ν

2

∇ enN+1 2 + B( enN+1 , u(tn+1 ),  enN+1 ) ≥ 0.

Proof. By using (3.5), we obtain

 n+1 1   n+1 ∇ u  ≤

fNj 2−1 . u0N 2 + Ca0 a0 ν j=1

Under the assumption of (3.11), we have

 n+1 1   ν

u0N 2 + Ca0 fNj 2−1 ≤ . a0 ν 2C0 j=1

That is

ν 2

− C0 ∇ u(tn+1 ) ≥ 0.

Then we obtain

ν

2

∇ enN+1 2 + B( enN+1 , u(tn+1 ),  enN+1 ) ≥

ν 2

− C0 ∇ u(tn+1 )

∇ enN+1 2 ≥ 0. 

Lemma 3.2.

a0 ≤ 2(1 − α )T α . bn

(3.12)

Proof.

bn = ((n + 1 )1−α − (n )1−α ) = n1−α

= n 1 −α

≥ n 1 −α

= n 1 −α

(1 − a ) n

(1 − a ) n

(1 − a ) 2n

+ + +



1+

(1 − α )(−α ) 1 n2

2!

(1 − α )(−α ) 1

2n

1−α

+

 −1



+ ···



n2

2!

(1 − a )

1 n

(1 − α )(−α ) 1



n2

2!

.

It is noticed that

(1 − a ) 2n

+

(1 − α )(−α ) 1 2!

n2

≥ 0.

Therefore, we have

2nα t α (2 − α ) a0 ≤ ≤ 2(1 − α )T α . bn (1 − α ) We show the error analysis for the solution of the scheme (3.3) and (3.4) in following theorem. Theorem 3.2. Let u be the solution of (1.1)–(1.3), Lemma (3.1), we have:

{unN }M−1 n=1



be the discrete solution of (3.3) and (3.4). Under the assumption of

√

u(·, tn ) − unN  + να∇ u(·, tn ) − ∇ unN  ≤ C (t 2−α + t −α N1−s ), n = 0, 1, . . . , M. Proof. For n = 0, the Eq. (3.3) can be written as

1 1 (u − u0N , ) + (∇ u1N , ∇ ) + B(u1N , u1N , ) + (∇ p1 , ) + (∇ u1N , q ) = f1N , . a0 N Subtracting (3.3) from a reformulation of (1.1) at t1 , we obtain

 ( e1N −  e0N , ) + a0 ν (∇  enN , ∇ ) + a0 B(u(t1 ), u(t1 ), ) − B(u1 , u1 , ) + a0 (∇ p(t1 ) − ∇ pN (t1 ), ) − a0 (∇ q, e1N ) = a0 (r 1 , ) + ((πN1 − I )(u(t1 ) − u(t0 )), ) + a0 ν ((πN1 − I )∇ u(t1 ), ∇ ).

(3.13)

486

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

Let = 2 e1N , q = 2( p(t1 ) − pN (t1 )), we have

2 e1N 2 + 2a0 ν∇  e1N 2 + B( e1N , u(t1 ),  e1N ) = a0 (r n+1 , 2 enN ) + 2((πN1 − I )(u(t1 ) − bn u(t0 )),  e1N ) + 2a0 (∇ p(t1 ) − ∇ p1N ,  e1N ) + 2a0 ν ((πN1 − I )∇ u(t1 ), ∇  e1N ) − a0 B(u1N ,  e1N ,  e1N ) − a0 B( e1N , u(t1 ),  e1N )





1 1 1 2  e1 2 + 2(πN1 − I )(u(t1 ) − u(t0 ))2 +  e  2a0 N 2 N  C ν ν a0 1 2 + a0 6(πN1 − I )∇ u(t1 )2 + ∇  e1N 2 + ∇ e1N 2 + ∇ eN  6 a0 6

≤ a0 2 a0 r 1 2 +



+ a0 C  e1N 2 +

ν

6

 ∇ e1N 2 + a0  p(t1 ) − p1N 2 + ∇  e1N 2

that is

 e1N 2 + a0 ν∇  e1N 2 ≤ C (a20 t 4−2α + Suppose

 ekN 2 + a0 ν∇  ekN 2 ≤ C

a20 t 4−2α b2k−1

1 2−2s N ). a0

+

N 2−2s bk−1 a0

 ,

k = 2, 3, . . . , n.

(3.14)

Next, we need to prove that it holds also for k = n + 1. Taking = 2 enN+1 and q = 2( p(tn+1 ) − pN (tn+1 )) in (3.3) and (3.4), we have:

2 enN+1 2 + 2a0 ν∇  enN+1 2 + B( enN+1 , u(tn+1 ),  enN+1 ) = 2a0 (r n+1 ,  enN+1 ) + 2((πN1 − I )(u(tn+1 ) −

n 

(b j − b j+1 )u(tn− j ) − bn u(t0 )),  enN+1 )

j=0

+ 2a0 ν ((πN1 − I )∇  enN+1 , ∇  enN+1 ) − 2a0 B(unN+1 ,  enN+1 ,  enN+1 ) − 2a0 B( enN+1 , u(tn+1 ),  enN+1 )



n−1 

+2 ≤ a0

a



(b j −

b j+1  enN− j

)

+

b n e0N ,  enN+1

j=0 0

bn

rn+1 2 +

+ 2a0 (∇ p(tn+1 ) − ∇ pnN+1 ,  enN+1 )



 bn 2  en+1 2 + (πN1 − I )(u(tn+1 ) − (b j − b j+1 )u(tn− j ) − bn u(t0 ))2 2a0 N bn n

j=0





bn n+1 2 ν +  e  + a0 8(πN1 − I )∇ u(tn+1 )2 + ∇  enN+1 2 + 2 N 8



+ a0 C  enN+1 2 + + a0



 p(tn+1 ) −

ν 8



∇ enN+1 2 +

n−1 

C a0

∇ enN+1 2 +

ν a0 8

∇ enN+1 2



(b j − b j+1 )( enN− j 2 +  enN+1 2 )

j=0

 + ∇

pnN+1 2

 .

 enN+1 2

Using the assumption (3.14) and the fact that bn−1− j > bn , we obtain

 enN+1 2 + a0 ν∇  enN+1 2 ≤C ≤C

=C =C

 a2 t 4−2α 0

bn

 a2 t 4−2α 0

bn

 a2 t 4−2α 0

b2n

 a2 t 4−2α 0

b2n

+

+

+

N 2−2s a0 N 2−2s a0 N 2−2s bn a0

N 2−2s + bn a0

Combing with (3.12), the estimate (3.13) is proved.



+C

n−1 

(b j − b j+1 )

 a2 t 4−2α 0

b2n−1− j

j=0



+C



n−1 

(b j − b j+1 )

 a2 t 4−2α

j=0

bn +

n−1 

0



b2n

+

N 2−2s bn−1− j a0

+

N 2−2s bn a0





(b j − b j+1 )

j=0



4. Numerical results In this part, we present numerical experiments to confirm the accuracy of the numerical solution obtained by the schemes (3.3) and (3.4). First, we need an exact solution to evaluate the accuracy of the numerical solution. Let ν = 1

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

487

Table 1 L2 -norm, H1 -norm, L∞ -norm error of u and p(α = 0.5).

u

p

t

L2 -norm

H1 -norm

L∞ -norm

0.1 0.01 0.001 0.0 0 01 0.1 0.01 0.001 0.0 0 01

9.7968E−003 7.3825E−004 3.5707E−005 1.3365E−006 0.5627 1.1980E−002 4.9805E−004 2.3293E−005

8.1088E−002 5.6148E−003 2.5551E−004 9.8406E−006 4.1180 0.1079 1.2636E−002 1.0175E−003

2.0903E−002 1.3658E−003 5.2698E−005 1.4110E−006 0.9819 7.8537E−002 1.1652E−002 9.5515E−004

Order 1.1229 1.3154 1.4268 1.6718 1.3812 1.3300

Table 2 L2 -norm, H1 -norm, L∞ -norm error of u and p(α = 0.9).

u

p

t

L2 -norm

H1 -norm

L∞ -norm

0.1 0.01 0.001 0.0 0 01 0.1 0.01 0.001 0.0 0 01

8.6206E−003 3.7258E−004 2.4443E−005 1.9278E−006 0.5459 6.4445E−003 9.6387E−005 8.3196E−005

6.9014E−002 2.6185E−003 1.7043E−004 1.3387E−005 4.0230 0.1081 5.3671E−002 2.0462E−004

1.5875E−002 4.4880E−004 2.8518E−005 2.2161E−006 0.9515 9.7875E−002 4.8977E−003 1.6358E−004

Order 1.3643 1.1831 1.1031 1.9280 1.8252 1.0639

Fig. 1. Spatial L2 , H1 and L∞ errors at α = 0.5, t = 0.0 0 01.

and the exact solution (u, p) of (1.1) and (1.2) to be

u(x, y, t ) =

 π t 2 sin2 (2π x ) sin(2π y ) − sin(2π x ) sin2 (π y ) ,

p(x, y, t ) = t 2 cos(π x ) sin(π y ). Then the associated forcing term f is given by f = Dtα u − u + u · ∇ u + ∇ p. Let unN , pnN are represented at the Lagrangian polynomials based on the Legendre–Gauss–Lobatto points ξ i , Legendre–Gauss nodes ηi , respectively,

unN (x, y ) =

N−1  N−1 

uni, j hi (x )h j (y ),

i=1 j=1

pnN (x, y ) =

N  N  i=0 j=0

pni, j li (x )l j (y ),

488

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

Fig. 2. First-component and second-component of the velocity errors at α = 0.5, t = 0.001.

Fig. 3. First-component and second-component of the velocity errors at α = 0.9, t = 0.001.

Fig. 4. The pressure errors at α = 0.5(left ), α = 0.9(right ), t = 0.001.

where uni, j = unN (ξi , ξ j ), pni, j = pnN (ηi , η j ) are the unknowns coefficient for discrete solution, li (x), hi (x) ∈ SN denote the Lagrange polynomial associated to Gauss–Lobatto points. The stop criteria of the algorithm until the following condition is satisfied

unN+1 − unN  ≤ 10−10 . The full discrete equations (3.3) and (3.4) are solved in  = (−1, 1 )2 with T = 1, N = 64. For 0 < α < 1 the convergence order of schemes (3.3) and (3.4) is O (t 2−α + t −α N 1−s ). In Tables 1 and 2, we compute the errors in three discrete norms and the temporal convergence rate for α = 0.5, 0.9, respectively. It shows that for α = 0.5, α = 0.9, the convergence orders of the velocity u are approximately 1.5 and 1.1, respectively, and at least first-order accurate for the pressure p. Because we need to solve a coupled nonlinear system for every time step, we cannot calculate for a sufficiently small amount of t. The L2 , H1 and L∞ errors at t = 0.0 0 01, α = 0.5 for u are shown in Fig. 1, it is seem that the error of u is influenced by the error in time direction. In Figs 2 and 3, we plot first-component and second-component of the velocity errors for the schemes schemes (3.3) and (3.4) at α = 0.5, 0.9. The above numerical tests agreement our theoretical analysis in Theorem 3.2. Although we do not get an estimate of the pressure, errors of the pressure is listed in Fig. 4. 5. Conclusion This article constructed a efficient numerical scheme to solve TFNSE. The numerical method is combined with a finite difference scheme in Caputo fractional derivative and Legendre-spectral method in space. We have derived a discrete stability inequality and error estimate for the numerical scheme which leads to 2 − α -order accuracy in time. Numerical examples are agreement with the theoretical prediction.

J. Zhang, J. Wang / Applied Mathematics and Computation 336 (2018) 481–489

489

Acknowledgment We are thankful to the reviewers for their useful comments which improved the quality of this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

C. Friedrich, Relaxation and retardation functions of the maxwell model with fractional derivatives, Rheol. Acta 30 (2) (1991) 151–158. W.G. Glöckle, T.F. Nonnenmacher, Fractional relaxation and the time-temperature superposition principle, Rheol. Acta 33 (4) (1994) 337–343. R. Hilfer, Applications of Fractional Calculus in Physics, 21, World Scientific, 20 0 0. D.Y. Song, T.Q. Jiang, Study on the constitutive equation with fractional derivative for the viscoelastic fluids-modified Jeffreys model and its application, Rheol. Acta 37 (5) (1998) 512–517. M. Li, J. Wang, Exploring delayed Mittag-Leffler type matrix functions to study finite time stability of fractional delay differential equations, Appl. Math. Comput. 324 (2018) 254–265. J.H. He, Approximate analytical solution for seepage flow with fractional derivatives in porous media, Comput. Methods Appl. Mech. Eng. 167 (1–2) (1998) 57–68. Z. Cao, J. Zhao, Z. Wang, F. Liu, L. Zheng, MHD flow and heat transfer of fractional maxwell viscoelastic nanofluid over a moving plate, J. Mol. Liq. 222 (2016) 1121–1127. J. Zhao, L. Zheng, X. Chen, X. Zhang, F. Liu, Unsteady Marangoni convection heat transfer of fractional maxwell fluid with Cattaneo heat flux, Appl. Math. Model. 44 (2017) 497–507. J. Zhao, L. Zheng, X. Zhang, F. Liu, Convection heat and mass transfer of fractional MHD maxwell fluid in a porous medium with Soret and Dufour effects, Int. J. Heat Mass Transf. 103 (2016) 203–210. J. Zhao, L. Zheng, X. Zhang, F. Liu, Unsteady natural convection boundary layer heat transfer of fractional maxwell viscoelastic fluid over a vertical plate, Int. J. Heat Mass Transf. 97 (2016) 760–766. J. Zhao, L. Zheng, X. Zhang, F. Liu, Unsteady natural convection heat transfer past a vertical flat plate embedded in a porous medium saturated with fractional Oldroyd-B fluid, J. Heat Transf. 139 (1) (2017). G.A. Baker, Galerkin Approximations for the Navier–Stokes Equations, Harvard University, August, 1976. D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2) (2001) 464–499. C. Min, F. Gibou, A second order accurate projection method for the incompressible Navier–Stokes equations on non-graded adaptive grids, J. Comput. Phys. 219 (2) (2006) 912–929. H. Okamoto, On the semi-discrete finite element approximation for the nonstationary Navier–Stokes equation, J. Fac. Sci. Univ. Tokyo Sect. A Math. 29 (3) (1982) 613–651. J. Shen, On error estimates of projection methods for Navier–Stokes equations: first-order schemes, SIAM J. Numer. Anal. 29 (1) (1992) 57–77. J. Shen, On error estimates of the projection methods for the Navier–Stokes equations: second-order schemes, Math. Comput. Am. Math. Soc. 65 (215) (1996) 1039–1065. F. Tone, Error analysis for a second order scheme for the Navier–Stokes equations, Appl. Numer. Math. 50 (1) (2004) 93–119. M. El-Shahed, A. Salem, On the Generalized Navier–Stokes Equations, Elsevier Science Inc., 2004. P.M.D. Carvalhoneto, G. Planas, Mild solutions to the time fractional Navier–Stokes equations in rn, J. Differ. Equ. 259 (7) (2015) 2948–2980. Y. Zhou, L. Peng, On the time-fractional Navier–Stokes equations, Comput. Math. Appl. 73 (2016) 874–891. Y. Zhou, L. Peng, Weak solutions of the time-fractional Navier–Stokes equations and optimal control, Comput. Math. Appl. 73 (6) (2017) 1016–1027. Z.Z. Ganji, D.D. Ganji, A.D. Ganji, M. Rostamian, Analytical solution of time-fractional Navier–Stokes equation in polar coordinate by homotopy perturbation method, Numer. Methods Part. Differ. Equ. 26 (1) (2010) 117–124. D. Kumar, J. Singh, S. Kumar, A fractional model of Navier–Stokes equation arising in unsteady flow of a viscous fluid, J. Assoc. Arab Univ. Basic Appl. Sci. 17 (2015) 14–19. S. Kumar, D. Kumar, S. Abbasbandy, M.M. Rashidi, Analytical solution of fractional Navier–Stokes equation by using modified laplace decomposition method, Ain Shams Eng. J. 5 (2) (2014) 569–574. X. Li, X. Yang, Y. Zhang, Error estimates of mixed finite element methods for time-fractional Navier–Stokes equations, J. Sci. Comput. 70 (2) (2017) 500–515. S. Momani, Z. Odibat, Analytical solution of a time-fractional Navier–Stokes equation by Adomian decomposition method, Appl. Math. Comput. 2 (177) (2006) 488–494. Y.M. Lin, C.J. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2) (2007) 1533–1552. X.J. Li, C.J. Xu, A space-time spectral method for the time fractional diffusion equation, Siam J. Numer. Anal. 47 (3) (2009) 2108–2131. A.H. Bhrawy, M.A. Zaky, D. Baleanu, New numerical approximations for space-time fractional Burgers’ equations via a Legendre spectral-collocation method, Roman. Rep. Phys. 67 (2015) 1–11. A.H. Bhrawy, M.A. Zaky, R.A.V. Gorder, A space-time Legendre spectral tau method for the two-sided space-time Caputo fractional diffusion-wave equation, Numer. Algorithms 71 (1) (2016) 151–180. A.H. Bhrawy, M.A. Zaky, J.T. Machado, Efficient Legendre spectral tau algorithm for solving two-sided space-time Caputo fractional advection-dispersion equation, J. Vib. Control 22 (8) (2016) 2053–2068. R. Temam, Navier–Stokes equations and nonlinear functional analysis, Soc. Ind. Appl. Math. 66 (1995). F. Liu, P. Zhuang, V. Anh, I. Turner, K. Burrage, Stability and convergence of the difference methods for the space time fractional advection diffusion equation, Appl. Math. Comput. 191 (1) (2007) 12–20. A.M. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, 1998. Y.M. Lin, X. Li, C.J. Xu, Finite difference/spectral approximations for the fractional cable equation, Math. Comput. 80 (275) (2011) 1369–1396.