Starting flows of a viscous incompressible fluid, modeled by the time-fractional derivatives, within a rotating channel due to an impulsive pressure gradient are studied. Using the eigenfunction expansion, the analytic solutions in series form are obtained. The flow of the ordinary fluid is studied as a special case of the time-fractional problem. The convergence of series solutions is proved. In addition, using the classical analytical method, coupled with the Laplace transform and Stehfest’s algorithm, an approximate solution is found. The flow rates in x- and y-directions are determined. In the case of the ordinary fluid, the steady-state and transient components of velocities are obtained. The numerical calculations are carried out by using the Mathcad software. It is found that, for fractional fluids, the reversal flow is much attenuated if the values of the fractional parameter are less than 1.

## I. INTRODUCTION

Fluid flows in rotating frames are important in many industrial applications. It is known that some natural phenomena, such as ocean circulations, hurricanes, tornadoes, and geophysical systems, imply rotating flows with heat and mass transfer. Theoretical and practical applications of the flows in rotating frames have led the researchers to attempt various problems related to such flows.

The starting flow in a rotating parallel channel due to an impulsive pressure gradient was conducted by Wang,^{1} using the separation of variables method and the Laplace transformation approach. Considering the Laplace transformation method with the homotopy perturbation method, Imran *et al.*^{2} have studied the influence of the Ekman number on the flow of a viscous fluid near a vertical plate, which is oscillating vertically and has a constant temperature. The unsteady Couette flow was studied by Das *et al.*^{3} in a system that was rotating. Hayat *et al.*^{4} have studied the magneto-hydrodynamic (MHD) transient flow of Maxwell fluids in a porous medium and rotating frame. An interesting study regarding electro-osmotic flows in a channel constituted by two plates placed parallel to each other or over an infinite plate in a rotating frame was conducted by Chang *et al*.^{5} Krishna *et al.*^{6} and Alghamdi *et al.*^{7,8} have analyzed the MHD convection flow with Hall effects through porous media in a rotating channel in the existence of a heat source dependent on the temperature. Other interesting problems with different boundary conditions can be found in Refs. 9–15.

Fractional calculus represents the more general differential calculus, which employs the derivatives/integrals of non-integer order. Not long ago, the differential calculus of non-integer order was found to be more appropriate to modeling of many practical processes.^{16}

In many areas, such as bioengineering, viscoelasticity, and polymer chemistry, fractional calculus provides much adequate models to physical systems compared to the models described by the classical ordinary differential equations.

It is known that the memory and hereditary properties of the materials are described by the time-fractional and space-fractional derivatives.^{17} Mathematical models illustrated by the fractional differential equations are used in electrical circuits, fractal theory, electromagnetic theory, etc.^{18} Some useful results can be found in Refs. 19 and 20.

The purpose of the present research article is to study a case not considered before, namely, the starting flow of a viscous incompressible fluid, with time-fractional derivatives, within a rotating channel due to an impulsive pressure gradient. Mainly, we want to study the impact of the fractional parameter on the flow parameters compared with the fluid modeled by the classical derivative (the case containing the fractional parameter α = 1).

The model characterized by the fractional differential equations is obtained from the model with the ordinary differential equations, by making use of the Caputo time-fractional derivatives instead of the classical time-derivatives.^{21} Using the eigenfunction expansion of the velocity field, the closed forms of the velocity components are obtained. These solutions are described by the Mittag-Leffler functions of two parameters.^{22} The convergence of the series solutions is proved. The case corresponding to the ordinary fluid (the fractional parameter α = 1) is studied as a particular case, and the steady-state and transient components of the velocity field are obtained. The Laplace transform of the complex velocity field is obtained by means of the classical analytical method in order to analyze the numerical results; however, although the inverse Laplace transform may be obtained, its mathematical expression is quite complicated and makes it difficult to use for numerical calculations. For this reason, an approximate solution of the velocity field is obtained using Stehfest’s algorithm for the inverse Laplace transform. The flow rates in the x-direction and y-direction are obtained. Graphs are used to illustrate the physical aspects of the fluid flow.

## II. PROBLEM STATEMENT

Consider the viscous fluid flow through a rectangular channel bounded by two rigid plates situated into planes z = 0 and z = h > 0, respectively, of a Cartesian coordinate system O_{xyz} (Fig. 1).

The system channel-fluid rotates with constant angular velocity, $\Omega \u2192=\Omega k\u2192$, *Ω* > 0, around the fixed z axis. In the present paper, we are looking for a velocity field of the form $v\u2192(z,t)=(u(z,t),v(z,t),0)$. It is found that the velocity field satisfies the continuity equation. The fluid motion is generated by a constant pressure gradient in the x-direction and by the channel rotation around the z axis. The time-fractional differential equations for a fluid flow in a rotating system are^{1,2,4,21}

where *λ* = *Ω*^{−1}, *ρ* is the density of the fluid, ν is the kinematic viscosity, G is the pressure gradient constant applied in the x-direction, H(t) is the unit step function, and $Dt\alpha c$ is the Caputo time-derivative operator of order *α* defined as^{17,18}

The appropriate initial and boundary conditions are

By adopting the following non-dimensional variables

we obtain the dimensionless problem (the star notation was dropped),

## III. SOLUTION OF THE PROBLEM

### A. Eigenfunction expansion of the velocity field

In this section, using the eigenfunction expansion,^{25–27} the solution of problems (6)–(8) is determined. In addition, some properties of the time-dependent velocity components are studied, and the convergence of the series solution is presented. Let *λ*_{n} = *n*^{2}*π*^{2}, $n\u2208N$, and $\psi n=2\u2009sin(n\pi z)$ the Dirichlet eigenvalues and eigenfunctions of the operator (−Δ) on the domain *D* = [0, 1] × (0, ∞).

The inner product $\u27e8f,\psi n\u27e9=\u222b01f(z)\psi n(z)dz$ in *L*^{2}(*D*) is used. Using eigenfunction decomposition, the solution *φ*(*z*, *t*) has the form

where the functions *φ*_{n}(*t*) satisfy the following conditions:

By using the Laplace transform to Eq. (10), with initial condition (11), and applying formula for the Laplace transform of the Caputo derivative,^{16,18} we obtain

where $\phi \xafn(s)=\u222b0\u221e\phi n(t)e\u2212stdt$ denotes the Laplace transform of function *φ*_{n}(*t*) and $an=n2\pi 2K2+2i$. Making use of the Mittag-Leffler function in two parameters,^{22,25} we have

respectively,

Some important properties of the function *φ*_{n}(*t*) are highlighted.

Functions

*φ*_{n}(*t*), n = 1, 2, … are continuous for*t*≥ 0 and*φ*_{n}(0) = 0.|

*φ*_{n}(*t*)| ≤*C*,*t*≥ 0,*C*=*const*.

We denote by *z*_{n} = −*t*^{α}*a*_{2n+1}. Using $an=n2\pi 2K2+2i$, we find $zn=\u2212t\alpha 2n+12\pi 2K2+2i$ and $\theta n=arg(zn)=\pi \u2212a\u2009tan2K2(2n+1)2\pi 2$, respectively.

It is apparent that *θ*_{n} ∈ (*π*/2, *π*) and *θ*_{n} < *θ*_{n+1}, n = 1, 2, ….

If $a\u2009tan2K29\pi 2<(2\u2212\alpha )\pi 2$, then $\alpha \pi 2<\theta 1$; therefore, there exists a real number μ, such that $\alpha \pi 2<\mu <\theta 1<\theta 2<\cdots <\pi .$ In the above conditions [see Ref. 23, Eq. (13) and Ref. 24, Eq. (2.10)],

where *p* ≥ 1 is an arbitrary integer and

For *t* → ∞, we have

and *φ*_{2n} = 0, and therefore, the functions *φ*_{n}(*t*) are bounded functions.

- (c)
The series solution given by Eq. (15) is uniformly convergent in D.

First, we note that functions {*ψ*_{n}(*z*)} are bounded functions on *z* ∈ [0, 1]. This property, together with property (b), implies that series (15) is uniformly convergent on D. Since all the terms used in the series are continuous functions on D, the sum *φ*(*z*, *t*) is a continuous function on D.

Equation (15) gives the closed form of the complex velocity field. To obtain the velocity components *u*(*z*, *t*) and *v*(*z*, *t*), the real and imaginary parts of *φ*(*z*, *t*) must be obtained. Denoting by $bn=(2n+1)2\pi 2K2,n=1,2,\u2026,K2=h2\Omega \nu $, it can be written

where $\alpha n=a\u2009tan2K2(2n+1)2\pi 2,\u2009n=1,2,\u2026$.

### B. Particular case α = 1(flow of the ordinary fluid)

In this case, by using the formula [see Ref. 25, Eq. (2.6)],

we obtain

The velocity components are given by

and can be written as the sum between the steady-state solution and the transient solution, as follows:

and

respectively. Our solutions, given by Eqs. (22)–(27), are equivalent to those obtained by Wang,^{1} if the rotation parameter *K*^{2} is replaced by 4β^{2}. In addition, it is worth highlighting that the steady-state solutions can be written in more suitable forms, namely,

where *d* = *K*/2.

### C. Flow rates

The volumetric flow rate Q is the volume of fluid that passes area *A* in the unit time. Let $n\u2192$ be the unit normal vector of area *A*, *da* the elementary area, and $u\u2192$ the velocity of fluid. The flow rate is the velocity integrated on the area, namely, $Q=\u222cAu\u2192\u22c5n\u2192da$.

For the flow in the x-direction, $u\u2192=uz,ti\u2192,\u2009n\u2192=i\u2192,\u2009da=dxdz,\u2009x\u2208[0,1],\u2009z\u2208[0,1]$, the flow rate is given by $Qx=\u222cAuz,tdxdz=\u222b0L\u222b01uz,tdxdz=L\u222b01uz,tdz$. Considering *L* = 1, the flow rate in the x-direction is given by

Using the velocity component given by Eq. (18), the rate of flow in the x-direction is

In a similar manner, we obtain the rate of flow in the y-direction given by

which, for the ordinary fluid, (*α* = 1), becomes

### D. Another form of the complex velocity. Numerical approximation

In this section, another form of the complex velocity is obtained by the classical analytic method together with the Laplace transformation method. Using the Laplace transform to Eqs. (6) and (8) with the initial condition (7), we obtain the following equations along with boundary conditions,

whose solution is

It is possible to obtain the inverse Laplace transform of function $\phi \xaf(z,s)$ given by Eq. (36) but the closed form of the corresponding complex velocity field $\phi (z,t)=L\u22121\phi \xaf(z,s)$ is rather complicated to be useful in numerical calculations. However, in the Appendix, we have presented the closed form of the inverse Laplace transform of the function $\phi \xaf(z,s)$.

There are many algorithms for the numerical calculation of the inverse Laplace transform. Stehfest’s formula, which approximates the inverse Laplace transform, is simple, and easy to use compared with other algorithms (for example, Talbot’s algorithm^{28}). Due to this reason, in the present paper, we use Stehfest’s algorithm^{29} for the inverse Laplace transform.

The approximate solution for the problem (6)–(8) is represented by

where k is a positive integer,

and $j+12$ is the integer part of the number (j + 1)/2.

## IV. NUMERICAL RESULTS AND DISCUSSION

To analyze the effect of the fractional order of the time-derivative on the velocity components *u*(*z*, *t*), *v*(*z*, *t*) and on the flow rates *Q*_{x}(*t*), *Q*_{y}(*t*), we performed some numerical calculations using the convergent series solutions (15) for the fractional parameter *α* ∈ (0, 1) and (20) and (21) for *α* = 1, respectively. Expression (15) is not very favorable for the numerical calculations. To obtain the solution (15) into a form favorable for numerical treatment, we use the method suggested by Gorenflo *et al.*,^{30} on the computation of the Mittag-Leffler functions. In agreement with the proposed algorithm^{28} [Eq. (19)], we shall use for the numerical calculations of the values of the Mittag-Leffler function *E*_{α,α+1}(*z*) the formula

For the velocity components *u* and *v*, given by the real part and the imaginary part of the complex velocity (15), respectively, we truncated the series at 500 terms for a very good approximation. The numerical values obtained by this way are compared with those given by Stehfest’s algorithm, namely, Eq. (37). The comparative results are shown in Table I for t = 2, K^{2} = 5, *α* = 0.5a, and p = 10, and a good agreement between the numerical values is found.

z . | U(z, t) Eq. (37) . | U(z, t) Eq. (15) . |
---|---|---|

0 | 0 | 0 |

0.1 | 0.024 591 5 | 0.024 590 3 |

0.2 | 0.041 297 5 | 0.041 297 6 |

0.3 | 0.051 900 6 | 0.051 900 6 |

0.4 | 0.057 721 | 0.057 721 4 |

0.5 | 0.059 570 6 | 0.059 570 8 |

0.6 | 0.057 722 4 | 0.057 721 4 |

0.7 | 0.051 902 | 0.051 900 6 |

0.8 | 0.041 298 6 | 0.041 297 6 |

0.9 | 0.041 298 6 | 0.041 297 6 |

1 | 0.024 591 5 | 0.024 590 3 |

z . | U(z, t) Eq. (37) . | U(z, t) Eq. (15) . |
---|---|---|

0 | 0 | 0 |

0.1 | 0.024 591 5 | 0.024 590 3 |

0.2 | 0.041 297 5 | 0.041 297 6 |

0.3 | 0.051 900 6 | 0.051 900 6 |

0.4 | 0.057 721 | 0.057 721 4 |

0.5 | 0.059 570 6 | 0.059 570 8 |

0.6 | 0.057 722 4 | 0.057 721 4 |

0.7 | 0.051 902 | 0.051 900 6 |

0.8 | 0.041 298 6 | 0.041 297 6 |

0.9 | 0.041 298 6 | 0.041 297 6 |

1 | 0.024 591 5 | 0.024 590 3 |

The primary starting velocity profiles in the x-direction are shown in Fig. 2, for K^{2} = 25, three values of time t and several values of the fractional parameter*α*. It is seen that the fluid modeled by the fractional differential equations has a different behavior than the ordinary fluid (*α* = 1). For large values of time t, the ordinary fluid has a larger flow reverse near the center of the channel. For the fluid modeled by the fractional time-derivatives, the reverse flow is more attenuated. Profiles of the secondary flows are plotted in Fig. 3. Notice that the transverse flow is toward the negative sense of the *y* axis.

Figures 4–7 show the primary flow rate in the x-direction and the secondary flow rate in the y-direction for low values of the rotation parameter (Figs. 4 and 5) and high values of the rotation parameter (Figs. 6 and 7), respectively. We have considered the case of the fractional fluid with *α* = 0.6 and the ordinary fluid for *α* = 1. We should highlight that, regarding the fractional time-derivative, the primary flow rate is higher if the fractional parameter value is lower. In addition, for the ordinary fluid, the primary flow rate is oscillating, and the reverse flow exists. For a fluid with time-fractional derivatives, the values of the primary flow rate tend fairly quickly to a constant value, and the reverse flow not appears or is very small. The secondary flow rate is not oscillating and becomes constant in a short time.

A comparison between the steady-state solutions *u*_{ss}(*z*) given by Eqs. (23) and (28) and v_{s}(*z*) given by Eqs. (26) and (29) is presented in Fig. 8 for two values of the rotation parameter *K*, respectively. It is observed from Fig. 8 that a very good agreement between velocity values is obtained.

The mathematical model of the studied problem clearly shows that the influence of the channel rotation on the fluid flow is characterized by the rotation parameter K. In the case of a fixed channel, the angular velocity is Ω = 0. The profiles of the primary velocity u(z, t) for the rotating channel and fixed channel, respectively, are plotted in Fig. 9 for two values of time t and for three values of the fractional parameter. It can be seen from this figure that the channel rotation leads to a slower flow in the x-direction compared with the flow when the rotation is blocked. This is due to the Coriolis force, which acts as a resistive force.

The aim of Fig. 10 is to emphasize the profiles of the flow rate *Q*_{x}(*t*) in the case of flows in a fixed channel (*K* = 0). It is observed from Fig. 10 that the flow rate in the x-direction increases with the fractional parameter and is higher than the flow rate in the case of flows in a rotating frame.

To have a validation of the above solutions, we obtained another form of the solution of Eq. (6) with the initial boundary conditions (7) and (8), namely,

## V. CONCLUDING REMARKS

The present work studies a new problem, namely, the unsteady Poiseuille flow in a rotating channel considering the governing equation as a fractional differential equation. The velocity components of primary and secondary flows, as well as the flow rates, are series solutions. The general forms of solutions are expressed with the Mittag-Leffler functions, and for the numerical computation of the values of these functions, the method proposed by Dingfelder and Weideman^{28} is used. In addition, by combining the Laplace transform method with the classical analytic method, the exact solution in the transformed domain is obtained.

A secondary set of numerical results is obtained by Stehfest’s algorithm. The matching between the numerical results, obtained by the two methods, highlights the correctness of calculations. Studying the effect of the fractional time-derivatives on the fluid velocity, it is found that if the order of the derivative is less than 1, then the fluid flows faster, and the flow reversal is more attenuated than for the ordinary fluid.

The parameter *α*, which represents the order of fractional derivative, remains as a control parameter. This parameter may be chosen so that the theoretical results of the fractional model are close to the experimental results. Thus, for some practical problems, a suitable fractional model may be chosen. This convenient choice could lead to a good agreement between the numerical data provided by the theoretical model and by experiments.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

The authors declare that they have no competing interest.

### APPENDIX: THE SOLUTION OF $\phi $(*z*, *t*)

Equation (36) can be written as

where

and

The inverse Laplace transform of function $\phi \xaf1(s)$ is

where $E\alpha ,\beta (z)=\u2211n=0\u221ezn\Gamma (\alpha n+\beta ),\u2009L\u22121s\alpha \u2212\beta s\alpha \u2212c=t\beta \u22121E\alpha ,\beta (ct\alpha )$ denotes the two-parameter Mittag-Leffler function.

To determine the inverse Laplace transform of function $\phi \xaf2(z,s)$, we consider the auxiliary function

whose inverse Laplace transform is

Note that $\phi \xaf2(z,s)=\varphi \xaf(z,w(s))$, therefore, the inverse Laplace transform of $\phi \xaf2(z,s)$ is given by

where

In the above relation, *Θ*(*β*, −*α*, −*K*^{2}*xs*^{−α}) is the Wright function given by

Finally, the inverse Laplace transform of function $\phi \xaf(z,s)$ is given by

## REFERENCES

_{α,β}(z) and its derivative