Numerical Solutions of Unsteady Boundary Layer Flow with a Time-Space Fractional Constitutive Relationship

Abstract

In this paper, we develop a new time-space fractional constitution relation to study the unsteady boundary layer flow over a stretching sheet. For the convenience of calculation, the boundary layer flow is simulated as a symmetrical rectangular area. The implicit difference method combined with an L1-algorithm and shift Grünwald scheme is used to obtain the numerical solutions of the fractional governing equation. The validity and solvability of the present numerical method are analyzed systematically. The numerical results show that the thickness of the velocity boundary layer increases with an increase in the space fractional parameter $\gamma$. For a different stress fractional parameter $\alpha$, the viscoelastic fluid will exhibit viscous or elastic behavior, respectively. Furthermore, the numerical method in this study is validated and can be extended to other time-space fractional boundary layer models.

1. Introduction

In recent years, the application of fractional operators in viscoelastic fluids has become more widespread. The Navier-Stokes (N-S) equations, motion equations describing the conservation of momentum of viscous incompressible fluids, were originally proposed by Navier in 1827. Due to the complex flow behavior of viscoelastic fluids, the relationship between strain and stress can no longer be linear. Combining the non-local properties and memory properties of the fractional derivative, it is applied to describe the flow behavior of viscoelastic fluid. Caputo and Mainardi [1,2] found that the exact results agree well with the experimental results when the fractional operator described a viscoelastic material. Scholars later conducted more in-depth studies based on this finding. To reveal the flow behavior of viscoelastic fluids, several models have been proposed, such as the Oldroyd-B and Maxwell models. Friendrich [3] introduced the Riemann–Liouville fractional operators to establish the Maxwell double fractional model with relaxation parameters. Makris et al. [4] predicted the dynamic response of viscoelastic fluids by using a general boundary-element formulation. Palade et al. [5] proposed the modified fractional Maxwell model with linear viscoelasticity. In their experiment, Li et al. analyzed the characteristics of xanthan gum and tianjing gum with a fractional order model [6,7,8]. Tan and Xu studied a plane surface suddenly set into motion in a fractional viscoelastic fluid, and obtained an exact solution of a viscoelastic fluid model between parallel plates [9,10]. Zhu et al. [11] investigated the shear stress and velocity distribution characteristics of Maxwell fluids. The peristaltic flow of a fractional Maxwell fluid through a medium was examined by Tripathi et al. [12]. Vieru et al. [13] discussed the flow of the fractional Maxwell fluids between two walls perpendicular to a plate. Applying the exponential rational function method, the solutions of nonlinear partial differential equations in different physical problems were given by Bekir and Kaplan [14]. More analyses on the numerical method of fractional order equations were given in the literature [15,16,17].

For the models discussed above, the space fractional derivative or the time fractional derivative the fractional parameters were usually used respectively. Recently, the time-space fractional models have been applied as a more effective tool than single-time or space fractional models [18,19,20]. Chen et al. [21] investigated the unsteady boundary layer flow with a time-space fractional derivative on a moving surface. A novel time-space heat conduction model with Cattaneo–Christov upper-convective derivative flux was proposed by Liu et al. [22]. Povstenko [23] studied a space–time fractional heat conduction equation and, furthermore, analyzed the significance of the time-space fractional derivatives. Atanacković developed a new generalization on heat conduction problem by using the space–time fractional derivative [24]. Chen et al. [25] studied the fractional constitutive equations characterizing magnetohydrodynamics (MHD) flow and heat transfer in viscoelastic fluids. Jannelli et al. [26] obtained numerical and analytical solutions of the time-space fractional diffusion equation. In recent articles, Pan et al. studied the fractional boundary layer problems of the Blasius flow of viscoelastic fluid over a plate [27]. Zhang et al. [28] considered the MHD flow of a fractional Maxwell fluid driven by a variable pressure gradient and obtained analytical and numerical solutions. Yang and Zhu [29] exhibited the pipe flow of a viscoelastic fluid by introducing the double fractional Maxwell model. To date, as we all know, the time-space fractional model on the boundary layer problem is still sparse, and the theoretical analysis of the numerical methods is a particularly challenging subject. To the best knowledge of the authors, the boundary layer flow with the time-space fractional constitutive relationship has never been mentioned.

Based on the above discussion, this paper proposed a new space-time fractional constitutive relationship to study the unsteady boundary layer flow of viscoelastic fluids. The finite difference method with the shifted Grünwald formulae and the L1 algorithm was used to obtain the numerical solution of the fractional governing equations. In addition, we also discussed the unique solvability of discrete equations and the validity of numerical methods. Finally, we analyzed the influence of all the parameters on the velocity distribution.

2. Governing Equations

The flow and heat transfer processes in the boundary layer are complex. As far as we know, the time-space fractional constitutive relationship on the boundary layer problem is sparse. Motivated by the above-mentioned discussions [3,16,17,21], the generalized law of Newton inner friction for the viscoelastic fluid used in this study can be expressed as:

$${\sigma }_{xy}+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }{\sigma }_{xy}}{\partial t^{\alpha }}={\mu }_{\gamma }\frac{{\partial }^{\gamma }u}{\partial y^{\gamma }},0\le \alpha \le 1,0\le \gamma \le 1,$$

(1)

where $\lambda$ is the relaxation time, ${\sigma }_{xy}$ is the shear stress, $\alpha$ and $\gamma$ are the fractional parameters, and ${\mu }_{\gamma }$ is the fractional dynamic viscosity, which depends on the properties of the fluid. ${\partial }^{\alpha }/\partial t^{\alpha }$ and ${\partial }^{\gamma }/\partial y^{\gamma }$ are the Caputo fractional order operator, and R–L are the fractional order operator [30], respectively, and are defined as:

$$\frac{{\partial }^{\alpha }{\sigma }_{xy}}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\xi \right)}^{-\alpha }\frac{\partial {\sigma }_{xy}\left(x,y,\xi \right)}{\partial \xi }}d\xi ,$$

(2)

$$\frac{{\partial }^{\gamma }u}{\partial y^{\gamma }}=\frac{1}{\Gamma \left(1-\gamma \right)}\frac{\partial }{\partial y}{{\displaystyle {\int }_o^y(y-\eta )}}^{-\gamma }u(x,\eta ,t)d\eta ,$$

(3)

where $\Gamma (·)$ is the Gamma function.

The two-dimensional unsteady boundary layer flow of viscoelastic fluids are considered over a stretching sheet with $u_w=a{x}^n$, where $a$ and $x$ are constants, and the coordinates are aligned along the stretching surface shown in Figure 1. $u$ and $v$ are the velocity components of the fluid in the $x$ and $y$ directions. At time $t=0$, the fluid and sheet are at rest, $u=0$ and $v=0$. When $t>0$, the fluid starts to flow due to the continuous stretching of the sheet with a linear velocity of $u_w=a{x}^n$. By introducing the constitutive Equation (1) to the basic equations of mass and momentum, we obtain the fractional boundary layer governing equations:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$$

(4)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+{\lambda }^{\alpha }\frac{{\partial }^{\alpha +1}u}{\partial t^{\alpha +1}}+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(u\frac{\partial u}{\partial x}\right)+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(v\frac{\partial u}{\partial y}\right)=\frac{1}{\rho }\frac{\partial }{\partial y}\left({\mu }_{\gamma }\frac{{\partial }^{\gamma }u}{\partial y^{\gamma }}\right)$$

(5)

For the convenience of the calculation, we introduce the following dimensionless variables:

$$u^{\ast }=\frac{u}{u_{\infty }},v^{\ast }=\frac{v}{u_{\infty }},x^{\ast }=\frac{x}{L_{\delta }},y^{\ast }=\frac{y}{L_{\delta }},t^{\ast }=\frac{t{u}_{\infty }}{L_{\delta }},{\lambda }^{\ast }=\frac{\lambda u_{\infty }}{L_{\delta }},\mathrm{and}\tilde{\mathrm{R}}\mathrm{e}=\frac{L_{\delta }{}^{\gamma }{u}_{\infty }\rho }{{\mu }_{\gamma }}.$$

The boundary layer governing equations become as:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$$

(6)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+{\lambda }^{\alpha }\frac{{\partial }^{1+\alpha }u}{\partial t^{1+\alpha }}+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}(u\frac{\partial u}{\partial x})+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}(v\frac{\partial u}{\partial y})=\frac{1}{\tilde{\mathrm{R}}\mathrm{e}}\frac{{\partial }^{1+\gamma }u}{\partial y^{1+\gamma }}$$

(7)

The corresponding initial and boundary conditions are:

$$t=0:u=0,v=0;$$

(8)

$$t>0:u=0,v=0\mathrm{at}x=0;u=a{x}^n,v=0\mathrm{at}y=0;u=0\mathrm{as}y\to \infty .$$

(9)

3. The Numerical Technique

Since the governing Equations (6) and (7) are highly coupled and nonlinear, to the best of our information, it is difficult to find their analytical solutions based on the current theory. Some numerical methods, and the calculation software have been applied as an effective tool to simulate the flow of non-Newtonian fluids [31,32,33]. In this section, the numerical solutions of governing Equations (6) and (7) are obtained by an effective finite difference method.

3.1. Discretization Method

In this section, the finite difference combined with the Grünwald scheme and L1-algorithm is employed to solve the governing equations. We defined $\Delta x$ and $\Delta y$ as special steps, and $\Delta t$ is the time step. $u_{i,j}^k$ and $v_{i,j}^k$ represent the values of the velocity of the fluid at the mesh point $\left(x_i,y_j,t_k\right)$.

First, for the Caputo fractional derivative ($0<\alpha \le 1$), we employ L1-approximation as follows [34]:

$$\frac{{\partial }^{\alpha }u(i,j,k)}{\partial t^{\alpha }}=\frac{\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[u(i,j,k)-{\displaystyle \sum _{s=1}^{k-1}\left({\alpha }_{s-1}^1-{\alpha }_s^1\right)u(i,j,k-s)-{\alpha }_{k-1}^{1}u(i,j,0)}\right]+O\left(\Delta t^{2-\alpha }\right)$$

(10)

For the R-L fractional derivative ($0<\gamma \le 1$), we employ the shifted Grünwald formulae as follows:

$$\frac{{\partial }^{\gamma +1}u\left(x_i,y_j,t_k\right)}{\partial y^{\gamma +1}}=\frac{1}{\Delta y^{\gamma +1}}{\displaystyle \sum _{l=0}^{j+1}\left(1-\frac{\gamma +2}{l}\right){\omega }_{l-1}{}^{\gamma +1}u\left(x_i,y_{j-l+1},t_k\right)}+O\left(\Delta y\right)$$

(11)

Second, we use the following discrete format for integer order terms:

$${\frac{\partial u}{\partial t}|}_{t=t_k}=\frac{u\left(x_i,y_j,t_k\right)-u\left(x_i,y_j,t_{k-1}\right)}{\Delta t}+O\left(\Delta t\right),$$

(12)

$${u\frac{\partial u}{\partial x}|}_{t=t_k}=u\left(x_i,y_j,t_{k-1}\right)\frac{u\left(x_i,y_j,t_k\right)-u\left(x_{i-1},y_j,t_k\right)}{\Delta y}+O\left(\Delta x+\Delta t\right),$$

(13)

$${v\frac{\partial u}{\partial y}|}_{t=t_k}=v\left(x_i,y_j,t_{k-1}\right)\frac{u\left(x_i,y_j,t_k\right)-u\left(x_i,y_{j-1},t_k\right)}{\Delta y}+O\left(\Delta y+\Delta t\right).$$

(14)

Applying the discrete format above, we have:

$${\frac{{\partial }^{\alpha +1}u}{\partial t^{\alpha +1}}|}_{t=t_k}=\frac{\Delta t^{-1-\alpha }}{\Gamma \left(2-\alpha \right)}\left[u_{i,j}^k-u_{i,j}^{k-1}-{\displaystyle \sum _{s=1}^{k-1}\left({\alpha }_{s-1}-{\alpha }_s\right)\left(u_{i,j}^{k-s}-u_{i,j}^{k-s-1}\right)}\right]+O(\Delta t),$$

(15)

$$\begin{array}{l}{\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(u\frac{\partial u}{\partial x}\right)|}_{t=t_k}=\frac{\Delta t^{-\alpha }}{\Delta x\Gamma \left(2-\alpha \right)}\left[u_{i,j}^{k-1}\left(u_{i,j}^k-u_{i-1,j}^k\right)-{\displaystyle \sum _{s=1}^{k-1}\left({\alpha }_{s-1}-{\alpha }_s\right)u_{i,j}^{k-s-1}\times \left(u_{i,j}^{k-s}-u_{i-1,j}^{k-s}\right)}\right]\\ +O\left(\Delta t+\Delta t^{2-\alpha }+\Delta x\right),\end{array}$$

(16)

$$\begin{array}{l}{\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(v\frac{\partial u}{\partial y}\right)|}_{t=t_k}=\frac{\Delta t^{-\alpha }}{\Delta y\Gamma \left(2-\alpha \right)}\left[v_{i,j}^{k-1}\left(u_{i,j}^k-u_{i,j-1}^k\right)-{\displaystyle \sum _{s=1}^{k-1}\left({\alpha }_{s-1}-{\alpha }_s\right)v_{i,j}^{k-s-1}\times \left(u_{i,j}^{k-s}-u_{i,j-1}^{k-s}\right)}\right]\\ +O\left(\Delta t+\Delta t^{2-\alpha }+\Delta y\right).\end{array}$$

(17)

Note that:

$$\begin{gathered} r_1={\lambda }^{\alpha }\frac{\Delta t^{-\alpha }}{\Gamma [2-\alpha ]};r_2=\frac{\Delta t}{\tilde{\mathrm{R}}\mathrm{e}\Delta y^{\gamma +1}};A_1={\displaystyle \sum _{s=1}^{k-1}[{\alpha }_{s-1}-{\alpha }_s](u_{i,j}^{k-s}-u_{i,j}^{k-s-1})}; \\ {A}_2={\displaystyle \sum _{s=1}^{k-1}[{\alpha }_{s-1}-{\alpha }_s]u_{i,j}^{k-s-1}(u_{i,j}^{k-s}-u_{i-1,j}^{k-s})};A_3={\displaystyle \sum _{s=1}^{k-1}[{\alpha }_{s-1}-{\alpha }_s]v_{i,j}^{k-s-1}(u_{i,j}^{k-s}-u_{i,j-1}^{k-s})}. \end{gathered}$$

Finally, the discrete format of the boundary layer governing equations is represented as follows:

$$v_{i,j}^k=v_{i,j-1}^k+\frac{\Delta y}{\Delta x}(u_{i-1,j}^k-u_{i,j}^k),$$

(18)

$$\begin{array}{l}-\left[\left(1+r_1\right)\frac{\Delta t}{\Delta y}{v}_{i,j}^{k-1}+r_2{\omega }_2^{\gamma +1}\right]u_{i,j-1}^k+\left[\left(1+r_1\right)+\frac{\Delta t}{\Delta x}\left(1+r_1\right)u_{i,j}^{k-1}+\left(1+r_1\right)\frac{\Delta t}{\Delta y}{v}_{i,j}^{k-1}-r_2{\omega }_1^{\gamma +1}\right]u_{i,j}^k\\ -r_2{\omega }_0^{\gamma +1}{u}_{i,j+1}^k-r_2{\displaystyle \sum _{l=0}^{j-2}{\omega }_{j-l+1}^{\gamma +1}{u}_{i,l}^k}=\left(1+r_1\right)u_{i,j}^{k-1}+r_1{A}_1+r_1\frac{\Delta t}{\Delta x}{A}_2+r_1\frac{\Delta t}{\Delta y}{A}_3+\left(1+r_1\right)\frac{\Delta t}{\Delta x}{u}_{i,j}^{k-1}{u}_{i-1,j}^k.\end{array}$$

(19)

The discrete formats of the initial and boundary conditions are written as:

$$u_{i,j}^0=v_{i,j}^0=0,u_{0,j}^k=v_{0,j}^k=0,u_{i,0}^k=a{x}^n,v_{i,0}^k=0,u_{i,N}^k=0.$$

3.2. Solvability of the Implicit Difference Scheme

In order to prove the unique solvability of the numerical equation, we rewrote Equations (19) and (20) in the matrix form as follows:

$$\begin{gathered} B{\mathit{V}}^k=-\frac{\Delta y}{\Delta x}B{\mathit{U}}^k, \\ G{U}^k=J{U}^{k-1}+\left(r_1\Delta t{\displaystyle \sum _{s=1}^{k-1}\left({\alpha }_{s-1}-{\alpha }_s\right)}\right)\times \left(I{\nabla }_t+{\Phi }_1{\nabla }_x+{\Phi }_2{\nabla }_y\right)U^{k-s} \end{gathered}$$

The coefficient matrix $G$ is:

$$G=\left[\begin{array}{cccccc}{H}_1& 0& 0& \cdots & 0& 0\\ Q_2& H_2& 0& \cdots & 0& 0\\ 0& Q_3& H_3& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & H_{M-1}& 0\\ 0& 0& 0& \cdots & Q_M& H_M\end{array}\right],\mathrm{B}=\left[\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\\ -1& 1& 0& \cdots & 0& 0\\ 0& -1& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& 0\\ 0& 0& 0& \cdots & -1& 1\end{array}\right],$$

where $H_i$ ($i=1,2,\cdots ,M$) and $Q_i$ ($i=2,\cdots ,M$) are the block matrixes:

$$\begin{gathered} H_i=\left[\begin{array}{cccccc}{a}_{i,1}& -r_2{\omega }_0^{\gamma +1}& 0& \cdots & 0& 0\\ b_{i,2}& a_{i,2}& -r_2{\omega }_0^{\gamma +1}& \cdots & 0& 0\\ -r_2{\omega }_3{}^{1+\gamma }& b_{i,3}& a_{i,3}& \cdots & 0& 0\\ -r_2{\omega }_4{}^{1+\gamma }& -r_2{\omega }_3{}^{1+\gamma }& b_{i,4}& \ddots & ⋮& ⋮\\ ⋮& ⋮& ⋮& \cdots & a_{i,N-2}& -r_2{\omega }_0^{\gamma +1}\\ -r_2{\omega }_{N-1}{}^{1+\gamma }& -r_2{\omega }_{N-2}{}^{1+\gamma }& -r_2{\omega }_{N-3}{}^{1+\gamma }& \cdots & b_{i,N-1}& a_{i,N-1}\end{array}\right],Q_i=-\frac{\Delta t}{\Delta x}(r_1+r_2)\left[\begin{array}{cccc}{u}_{i,1}^{k-1}& 0& \cdots & 0\\ 0& u_{i,2}^{k-1}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& u_{i,N-1}^{k-1}\end{array}\right], \\ \{\begin{array}{ll}{a}_{i,j}=\left(1+r_1\right)+\frac{\Delta t}{\Delta x}\left(1+r_1\right)u_{i,j}^{k-1}+\left(1+r_1\right)\frac{\Delta t}{\Delta y}{v}_{i,j}^{k-1}-r_2{\omega }_1^{\gamma +1}\hfill & ,j=1,2,\cdots ,N-1\\ b_{i,j}=-\left(1+r_1\right)\frac{\Delta t}{\Delta y}{v}_{i,j}^{k-1}-r_2{\omega }_2^{\gamma +1}\hfill & ,j=2,\cdots ,N-1\end{array} \end{gathered}$$

and the coefficient matrix $J$, ${\Phi }_1$, ${\Phi }_2$ are:

$$\begin{gathered} \mathrm{J}=\left[\begin{array}{cccc}1+r_1& 0& \cdots & 0\\ 0& 1+r_1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1+r_1\end{array}\right], \\ {\Phi }_1=diag\left(u_{1,1}^{k-s-1},u_{1,2}^{k-s-1},\cdots ,u_{1,N-1}^{k-s-1},u_{2,1}^{k-s-1},u_{2,2}^{k-s-1},\cdots ,u_{2,N-1}^{k-s-1},\cdots ,u_{M,1}^{k-s-1},u_{M,2}^{k-s-1},\cdots ,u_{M,N-1}^{k-s-1}\right), \\ {\Phi }_2=diag\left(v_{1,1}^{k-s-1},v_{1,2}^{k-s-1},\cdots ,v_{1,N-1}^{k-s-1},v_{2,1}^{k-s-1},v_{2,2}^{k-s-1},\cdots ,v_{2,N-1}^{k-s-1},\cdots ,v_{M,1}^{k-s-1},v_{M,2}^{k-s-1},\cdots ,v_{M,N-1}^{k-s-1}\right). \end{gathered}$$

where:

$$\begin{gathered} U^k={\left[u_{1,1}^k,u_{1,2}^k,\cdots u_{1,N-1}^k,u_{2,1}^k,u_{2,2}^k,\cdots ,u_{2,N-1}^k,\cdots ,u_{M,1}^k,\cdots u_{M,N-1}^k\right]}^{\mathrm{T}}, \\ {V}^k={\left[v_{1,1}^k,v_{1,2}^k,\cdots v_{1,N-1}^k,v_{2,1}^k,v_{2,2}^k,\cdots ,v_{2,N-1}^k,\cdots ,v_{M,1}^k,\cdots v_{M,N-1}^k\right]}^{\mathrm{T}}, \\ {U}^{k-1}={\left[u_{1,1}^{k-1},u_{1,2}^{k-1},\cdots u_{1,N-1}^{k-1},u_{2,1}^{k-1},u_{2,2}^{k-1},\cdots ,u_{2,N-1}^{k-1},\cdots ,u_{M,1}^{k-1},\cdots u_{M,N-1}^{k-1}\right]}^{\mathrm{T}}, \\ {U}^{k-s}={\left[u_{1,1}^{k-s},u_{1,2}^{k-s},\cdots u_{1,N-1}^{k-s},u_{2,1}^{k-s},u_{2,2}^{k-s},\cdots ,u_{2,N-1}^{k-s},\cdots ,u_{M,1}^{k-s},\cdots u_{M,N-1}^{k-s}\right]}^{\mathrm{T}}, \\ {\nabla }_t{u}_{i,j}^k=\frac{u_{i,j}^k-u_{i,j}^{k-1}}{\Delta t},{\nabla }_x{u}_{i,j}^k=\frac{u_{i,j}^k-u_{i-1,j}^k}{\Delta x},{\nabla }_y{u}_{i,j}^k=\frac{u_{i,j}^k-u_{i,j-1}^k}{\Delta y}. \end{gathered}$$

From the above, we can clearly observe the $G$ and $B$ are full rank matrices, then the numerical schemes Equations (19) and (20) are uniquely solvable.

4. Validation of the Numerical Discretization Method

To verify the correctness of the numerical discretization method, we construct new governing equations by introducing $f\left(x,y,t\right)$:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$$

(20)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+{\lambda }^{\alpha }\frac{{\partial }^{1+\alpha }u}{\partial t^{1+\alpha }}+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(u\frac{\partial u}{\partial x}\right)+{\lambda }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}\left(v\frac{\partial u}{\partial y}\right)=\frac{1}{\tilde{\mathrm{R}}\mathrm{e}}\frac{{\partial }^{1+\gamma }u}{\partial y^{1+\gamma }}+f\left(x,y,t\right).$$

(21)

with initial and boundary conditions, which are:

$$u\left(x,y,0\right)=v\left(x,y,0\right)=0,u\left(x,0,t\right)=v\left(x,0,t\right),u\left(x,1,t\right)=0$$

(22)

The exact solutions can be obtained as follows:

$$u\left(x,y,t\right)=x^2{\left(x-1\right)}^2{y}^2{\left(1-y\right)}^2{t}^2,$$

(23)

$$v\left(x,y,t\right)=2x\left(x-1\right)\left(1-2x\right)\left(\frac{1}{3}{y}^3-\frac{1}{2}{y}^4+\frac{1}{5}{y}^5\right)t^2.$$

(24)

Combining the initial and boundary conditions, we obtain the following expression of $f\left(x,y,t\right)$:

$$\begin{array}{l}f\left(x,y,t\right)=2x^2{\left(1-x\right)}^2{y}^2{\left(1-y\right)}^{2}t\left(1+\frac{{\lambda }^{\alpha }{t}^{-\alpha }}{\Gamma \left(2-\alpha \right)}\right)+2x^3{\left(1-x\right)}^3{y}^4{\left(1-y\right)}^4\left(1-2x\right)\left(t^4+\frac{24{\lambda }^{\alpha }{t}^{4-\alpha }}{\Gamma \left(5-\alpha \right)}\right)\\ -4x^3{\left(1-x\right)}^{3}y\left(1-y\right)\left(1-2x\right)\left(1-2y\right)\left(\frac{1}{3}{y}^3-\frac{1}{2}{y}^4+\frac{1}{5}{y}^5\right)\left(t^4+\frac{24{\lambda }^{\alpha }{t}^{4-\alpha }}{\Gamma \left(5-\alpha \right)}\right)\\ -\frac{2}{\tilde{\mathrm{R}}\mathrm{e}}{x}^2{\left(1-x\right)}^2{t}^2\left(\frac{y^{1-\gamma }}{\Gamma \left(2-\gamma \right)}+\frac{12y^{3-\gamma }}{\Gamma \left(4-\gamma \right)}-\frac{6y^{2-\gamma }}{\Gamma \left(3-\gamma \right)}\right).\end{array}$$

(25)

Equations (21) and (22) are solved numerically by applying the implicit difference method. Figure 2 shows a comparison of numerical and analytical solutions. As the figure shows, the two curves are almost identical, which can indicate the correctness of our numerical method. In order to test the convergence of the numerical method,

Table 1. Temporal error and convergence order for different $\gamma$ with $\alpha =0.1$.

$\alpha =0.1$ $\gamma =0.98$	$L_2$Error	Order	$L_{\infty }$Error	Order
1/8	3.9055e-05	-	8.4866e-05	-
1/16	1.9529e-05	0.9989	4.2644e-05	0.9282
1/32	9.7050e-06	1.0088	2.1337e-05	0.9901
1/64	4.8460e-06	1.0019	1.0655e-05	1.0018
$\alpha =0.1$ $\gamma =0.95$	$L_2$Error	Order	$L_{\infty }$Error	Order
1/8	4.0576e-05	-	8.9865e-05	-
1/16	2.0305e-05	0.9988	4.4954e-05	0.9993
1/32	1.0087e-05	1.0093	2.2335e-05	1.0092
1/64	5.0183e-06	1.0071	1.1127e-05	1.0052

and

Table 2. Temporal error and convergence order for different $\alpha$ with $\gamma =0.8$.

$\alpha =0.3$ $\gamma =0.8$	$L_2$Error	Order	$L_{\infty }$Error	Order
1/8	5.7726e-05	-	1.3641e-04	-
1/16	2.6869e-05	1.1033	6.3664e-05	1.0993
1/32	1.2973e-05	1.0504	3.0797e-05	1.0477
1/64	6.3905e-06	1.0215	1.5257e-05	1.0133
$\alpha =0.4$ $\gamma =0.8$	$L_2$Error	Order	$L_{\infty }$Error	Order
1/8	5.6212e-05	-	1.3151e-04	-
1/16	2.6040e-05	1.1101	6.1069e-05	1.1067
1/32	1.2539e-05	1.0543	2.9888e-05	1.0309
1/64	6.1704e-06	1.0230	1.4809e-05	1.0131

present the $L_2$ error, $L_{\infty }$ error, and convergence order of $\Delta y$ for different fractional parameters $\gamma$ and $\alpha$ at $t=1$. As indicated by the tabular data, the numerical results are in good agreement with the analytical solution, and the convergence order reaches the expected first order. Thus, we conclude that the numerical technique in this study is valid and can be extended to other time-space fractional boundary layer problems.

5. Results and Discussion

In the calculation, we choose $X_{\max }=1$ and $Y_{\max }=7$ ($X_{\max }$ represents a point on the surface that can be extended qualitatively, and $Y_{\max }$ corresponds to $y\to \infty$). The fractional governing equations at each internal node form a system of equations at each specific time node that can be solved by the Gauss algorithm. After multiple numerical simulations, we assign the space and time step sizes $\Delta x=0.05$, $\Delta y=0.05$, and $\Delta t=0.01$. Differently sized spatial grid steps were selected to test the convergence of the connection.

In this section, the effects of all the parameters embedded in the velocity and temperature within the boundary layer are graphically illustrated. Without loss of generality, the other parameters are fixed as $\alpha =0.1$, $\gamma =0.9$, $\lambda =0.1$, and $\tilde{\mathrm{R}}\mathrm{e}=5$ when examining the influence of a specific parameter. The results in Figure 3 show that the solutions for different grid steps are almost consistent, which indicates that the selected step size is valid. Figure 4 presents three-dimensional distributions of dimensionless velocity by one set of parameters, which demonstrate good stability and convergence of our numerical solutions.

Figure 5 and Figure 6 demonstrate the spatial evolution of the nondimensional velocity in the steady-state, where the value of $x$ is taken as $X_{\max }=1$ with space fractional parameter $\gamma$ and time fractional parameter $\alpha$. Figure 5 shows the distribution of the velocity over differences in the space fractional parameter $\gamma$ and indicates that the change of velocity versus the space fractional parameter $\gamma$ is consistent. The smaller the fractional parameter $\gamma$ is, the smaller the velocity is. That is, the thickness of the velocity boundary layer decreases with a decrease in space fractional parameter $\gamma$. The influence of time fractional parameter $\alpha$ on the velocity field is given in Figure 6. The smaller the value of $\alpha$, the faster the velocity decreases, accompanied by a decrease in the velocity boundary layer thickness.

Figure 7 shows the distribution of the velocity for a different stress fractional parameter $\alpha$ in response to $t$. As the stress fractional parameter $\alpha$ increases, it can be seen that the velocity will eventually appear as an oscillations phenomenon and then stabilize at a greater speed, i.e., the viscoelastic fluid presents anelasticity behavior. On the contrary, for a smaller stress fractional parameter $\alpha$, the fluid exhibits stronger viscous behavior, which stabilizes the flow in a faster time. Figure 8 presents the velocity distribution for different relaxation time parameter $\lambda$ in response to $t$. It can be found that as the relaxation parameter $\lambda$ increases, the velocity will stabilize in a longer time. In other words, the relaxation parameter results in a delay on the flow process. Therefore, the velocity tends to stabilize more slowly for a larger relaxation parameter $\lambda$.

6. Conclusions

This paper presents the flow of a fractional viscoelastic fluid over a stretching sheet. The space-time fractional constitutive relationship is introduced into the momentum equation. The shifted Grünwald formula is used to approximate the space fractional derivative, and an L1-scheme is applied to approximate the time-fractional derivative. Then the governing equations are solved numerically by an effective finite difference method. The effects of the space fractional parameter $\gamma$ and the time fractional parameter $\alpha$ in response to $y$ are analyzed. The numerical results indicate that the thickness of the velocity boundary layer increases with an increase in the space fractional parameter $\gamma$. Similarly, the smaller the fractional parameter $\alpha$ is, the thinner the velocity boundary layer is. In addition, we also analyzed the influence of stress fractional parameters $\alpha$ and the relaxation parameters $\lambda$ in response to $t$ on the velocity distribution. It is found that for a larger stress fractional parameter $\alpha$, the viscoelastic fluid exhibits a stronger elastic behavior, and for a smaller $\alpha$, the fluid exhibits a stronger viscous behavior. In addition, a larger relaxation parameter will cause a delay in the flow process.