Numerical simulations of unsteady flows in a rotating channel using a novel eigenfunction expansion method

Rubbab, Q.; Mahsud, Y.; Irshad, S.; Imran, M. A.; Ahmadian, A.; Salahshour, S.; Ferrara, M.

doi:10.1063/5.0012874

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,¹ using the separation of variables method and the Laplace transformation approach. Considering the Laplace transformation method with the homotopy perturbation method, Imran et al.² 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.³ in a system that was rotating. Hayat et al.⁴ 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.⁵ Krishna et al.⁶ 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.¹⁶

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.¹⁷ Mathematical models illustrated by the fractional differential equations are used in electrical circuits, fractal theory, electromagnetic theory, etc.¹⁸ 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.²¹ 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.²² 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).

FIG. 1.

View large Download slide

The flow geometry.

The system channel-fluid rotates with constant angular velocity, $\vec{Ω} = Ω \vec{k}$ ⁠, Ω > 0, around the fixed z axis. In the present paper, we are looking for a velocity field of the form $\vec{v} (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

λ^{α} {}_{c}D_{t}^{α} u - 2 v = \frac{λ}{ρ} G H (t) + ν λ \frac{\partial^{2} u}{\partial z^{2}},

(1)

λ^{α} {}_{c}D_{t}^{α} v + 2 u = ν λ \frac{\partial^{2} v}{\partial z^{2}},

(2)

where λ = Ω⁻¹, ρ 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 ${}_{c}D_{t}^{α}$ is the Caputo time-derivative operator of order α defined as^17,18

\begin{matrix} {}_{c}D_{t}^{α} f (z, t) = \{\begin{gathered} \int_{0}^{t} \frac{\partial f (z, τ)}{\partial τ} \frac{{(t - τ)}^{- α}}{Γ (1 - α)} d τ, α \in (0,1), \\ \frac{\partial f (z, t)}{\partial t}, α = 1 . \end{gathered}) \end{matrix}

(3)

The appropriate initial and boundary conditions are

\begin{matrix} \begin{gathered} u (z, 0) = v (z, 0) = 0, z \in [0, h], \\ u (0, t) = v (0, t) = 0, t \geq 0, \\ u (h, t) = v (h, t) = 0, t \geq 0 . \end{gathered} \end{matrix}

(4)

By adopting the following non-dimensional variables

t^{*} = \frac{t}{λ}, z^{*} = \frac{z}{h}, (u^{*}, v^{*}) = \frac{ν ρ (u, v)}{G h^{2}}, K^{2} = \frac{h^{2} Ω}{ν},

(5)

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

K^{2} {}_{c}D_{t}^{α} φ + 2 i K^{2} φ = H (t) + \frac{\partial^{2} φ}{\partial z^{2}},

(6)

φ (z, 0) = 0, z \in [0,1],

(7)

φ (0, t) = φ (1, t) = 0, t \geq 0,

(8)

where K² denotes the rotation parameter, φ(z, t) = u(z, t) + iv(z, t), and i² = −1 is the complex velocity field. Since the initial and boundary conditions (7) and (8) are homogeneous, Eq. (6) has a unique solution.

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²π², $n \in N$ ⁠, and $ψ_{n} = \sqrt{2} sin (n π z)$ the Dirichlet eigenvalues and eigenfunctions of the operator (−Δ) on the domain D = [0, 1] × (0, ∞).

The inner product $⟨ f, ψ_{n} ⟩ = \int_{0}^{1} f (z) ψ_{n} (z) d z$ in L²(D) is used. Using eigenfunction decomposition, the solution φ(z, t) has the form

φ (z, t) = \sum_{n = 1}^{\infty} φ_{n} (t) ψ_{n} (z),

(9)

where the functions φ_n(t) satisfy the following conditions:

K^{2} {}_{c}D_{t}^{α} φ_{n} (t) + 2 i K^{2} φ_{n} (t) = \frac{1 - {(- 1)}^{n}}{n π} \sqrt{2} H (t) - λ_{n} φ_{n} (t),

(10)

φ_{n} (0) = 0, n = 1,2, \dots .

(11)

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

{\bar{φ}}_{n} (s) = \sqrt{2} \frac{1 - {(- 1)}^{n}}{n π} \frac{s^{- 1}}{s^{α} + a_{n}},

(12)

where ${\bar{φ}}_{n} (s) = \int_{0}^{\infty} φ_{n} (t) e^{- s t} d t$ denotes the Laplace transform of function φ_n(t) and $a_{n} = \frac{n^{2} π^{2}}{K^{2}} + 2 i$ ⁠. Making use of the Mittag-Leffler function in two parameters,^22,25 we have

L^{- 1} \{\frac{s^{α - β}}{s^{α} + b}\} = t^{β - 1} E_{α, β} (- b t^{α}), E_{α, β} (z) = \sum_{k = 0}^{\infty} \frac{z^{k}}{Γ (α k + β)}, α, β > 0,

(13)

respectively,

φ_{n} (t) = \sqrt{2} \frac{1 - {(- 1)}^{n}}{n π K^{2}} t^{α} E_{α, α + 1} (- a_{n} t^{α}) .

(14)

Using (14) into Eq. (9), we obtain

\begin{matrix} φ (z, t) & = & \sum_{n = 1}^{\infty} \frac{\sqrt{2} [1 - {(- 1)}^{n}]}{n π K^{2}} t^{α} E_{α, α + 1} (- a_{n} t^{α}) ψ_{n} (z) \\ = & \frac{2 \sqrt{2}}{π K^{2}} \sum_{n = 0}^{\infty} \frac{t^{α} ψ_{2 n + 1} (z)}{2 n + 1} E_{α, α + 1} (- a_{2 n + 1} t^{α}) . \end{matrix}

(15)

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 $a_{n} = \frac{n^{2} π^{2}}{K^{2}} + 2 i$ ⁠, we find $z_{n} = - t^{α} [\frac{{(2 n + 1)}^{2} π^{2}}{K^{2}} + 2 i]$ and $θ_{n} = arg (z_{n}) = π - a tan (\frac{2 K^{2}}{{(2 n + 1)}^{2} π^{2}})$ ⁠, respectively.

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

If $a tan (\frac{2 K^{2}}{9 π^{2}}) < \frac{(2 - α) π}{2}$ ⁠, then $\frac{α π}{2} < θ_{1}$ ⁠; therefore, there exists a real number μ, such that $\frac{α π}{2} < μ < θ_{1} < θ_{2} < \dots < π .$ In the above conditions [see Ref. 23, Eq. (13) and Ref. 24, Eq. (2.10)],

t^{α} E_{α, α + 1} (z_{n}) = - \sum_{k = 1}^{p} \frac{t^{α}}{Γ (1 + α (1 - k)) z_{n}^{k}} + O ({|z_{n}|}^{- p - 1}),

(16)

where p ≥ 1 is an arbitrary integer and

\begin{matrix} lim_{t \to \infty} t^{α} E_{α, α + 1} (z_{n}) & = & \sum_{k = 1}^{p} \frac{{(- 1)}^{k + 1}}{Γ (1 + α (1 - k)) a_{2 n + 1}^{k}} lim_{t \to \infty} \frac{1}{t^{α (k - 1)}} \\ = & \frac{1}{a_{2 n + 1}} . \end{matrix}

For t → ∞, we have

|φ_{2 n + 1} (t)| \leq \frac{2 \sqrt{2}}{(2 n + 1) π K^{2}} \frac{1}{|a_{2 n + 1}|} < \frac{\sqrt{2}}{(2 n + 1) π K^{2}} < \frac{\sqrt{2}}{π K^{2}} vf,

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 $b_{n} = \frac{{(2 n + 1)}^{2} π^{2}}{K^{2}}, n = 1,2, \dots, K^{2} = \frac{h^{2} Ω}{ν}$ ⁠, it can be written

\begin{matrix} E_{α, α + 1} (- a_{2 n + 1} t^{α}) & = & E_{α, α + 1} [- (b_{n} + 2 i) t^{α}] = \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} t^{α k} {(b_{n} + 2 i)}^{k}}{Γ (1 + α (k + 1))} \\ = & \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} t^{α k} {(b_{n}^{2} + 4)}^{k / 2} [cos (k α_{n}) + i sin (k α_{n})]}{Γ (1 + α (k + 1))}, \end{matrix}

(17)

where $α_{n} = a tan (\frac{2 K^{2}}{{(2 n + 1)}^{2} π^{2}}), n = 1,2, \dots$ ⁠.

Using (17) into Eq. (15), we achieve the following expressions for the velocity components:

\begin{matrix} u (z, t) & = & R e [φ (z, t)] = \frac{2 \sqrt{2}}{K^{2}} \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} \frac{t^{α} ψ_{2 n + 1} (z)}{(2 n + 1) π} \\ \times \frac{{(- 1)}^{k} t^{α k} {(b_{n}^{2} + 4)}^{k / 2} cos (k α_{n})}{Γ [1 + α (k + 1)]}, \\ v (z, t) & = & I m [φ (z, t)] = \frac{2 \sqrt{2}}{K^{2}} \sum_{n = 0}^{\infty} \sum_{k = 0}^{\infty} \frac{t^{α} ψ_{2 n + 1} (z)}{(2 n + 1) π} \\ \times \frac{{(- 1)}^{k} t^{α k} {(b_{n}^{2} + 4)}^{k / 2} sin (k α_{n})}{Γ [1 + α (k + 1)]} . \end{matrix}

(18)

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

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

E_{1,2} (z) = \frac{1}{z} (e^{z} - 1),

(19)

we obtain

\begin{matrix} E_{1,2} (- t a_{2 n + 1}) & = & \frac{1}{t (b_{n} + 2 i)} [1 - e^{- t (b_{n} + 2 i)}] \\ = & \frac{[2 sin (2 t) - b_{n} cos (2 t)] e^{- b_{n} t} + b_{n}}{t (b_{n}^{2} + 4)} \\ + \frac{i \{[b_{n} s i n (2 t) + 2 c o s (2 t)] e^{- b_{n} t} - 2\}}{t (b_{n}^{2} + 4)} . \end{matrix}

The velocity components are given by

\begin{matrix} u (z, t) & = & \frac{2 \sqrt{2}}{π K^{2}} \sum_{n = 0}^{\infty} \frac{ψ_{2 n + 1} (z)}{(2 n + 1) (b_{n}^{2} + 4)} \\ \times \{[2 sin (2 t) - b_{n} cos (2 t)] e^{- b_{n} t} + b_{n}\}, \end{matrix}

(20)

\begin{matrix} v (z, t) & = & \frac{2 \sqrt{2}}{π K^{2}} \sum_{n = 0}^{\infty} \frac{ψ_{2 n + 1} (z)}{(2 n + 1) (b_{n}^{2} + 4)} \\ \times \{[b_{n} sin (2 t) + 2 cos (2 t)] e^{- b_{n} t} - 2\}, \end{matrix}

(21)

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

u (z, t) = u_{s s} (z) + u_{t} (z, t),

(22)

u_{s s} (z) = 2 \sqrt{2} π \sum_{n = 0}^{\infty} \frac{(2 n + 1) sin [(2 n + 1) π z]}{{(2 n + 1)}^{4} π^{4} + 4 K^{4}},

(23)

u_{t} (z, t) = \frac{2 \sqrt{2}}{π K^{2}} \sum_{n = 0}^{\infty} \frac{sin [(2 n + 1) π z]}{(2 n + 1) (b_{n}^{2} + 4)} [2 sin (2 t) - b_{n} cos (2 t)] e^{- b_{n} t},

(24)

and

v (z, t) = v_{s s} (z) + v_{t} (z, t),

(25)

v_{s s} (z) = \sum_{n = 0}^{\infty} \frac{- 4 K^{2} \sqrt{2} sin [(2 n + 1) π z]}{π (2 n + 1) [{(2 n + 1)}^{4} π^{4} + 4 K^{4}]},

(26)

v_{t} (z, t) = \frac{2 \sqrt{2}}{π K^{2}} \sum_{n = 0}^{\infty} \frac{sin [(2 n + 1) π z]}{(2 n + 1) (b_{n}^{2} + 4)} [b_{n} sin (2 t) + 2 cos (2 t)] e^{- b_{n} t},

(27)

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

\begin{matrix} u_{s s} (z) & = & \frac{\sqrt{2}}{4 K^{2}} \frac{cosh [(2 z - 1) d] sinh (d) sin (2 d z)}{cos h^{2} (d) - si n^{2} (d)} \\ - \frac{\sqrt{2}}{4 K^{2}} \frac{s i n [(2 z - 1) d] cos (d) sinh (2 d z)}{cos h^{2} (d) - si n^{2} (d)}, \end{matrix}

(28)

\begin{matrix} v_{s s} (z) & = & \frac{\sqrt{2}}{4 K^{2}} \frac{cosh [(2 z - 1) d] cosh (d) cos (2 d z)}{cos h^{2} (d) - si n^{2} (d)} \\ + \frac{\sqrt{2}}{4 K^{2}} (\frac{\sin [(2 z - 1) d] sin (d) cosh (2 d z)}{cos h^{2} (d) - si n^{2} (d)} - 1), \end{matrix}

(29)

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 $\vec{n}$ be the unit normal vector of area A, da the elementary area, and $\vec{u}$ the velocity of fluid. The flow rate is the velocity integrated on the area, namely, $Q = \iint_{A} \vec{u} \cdot \vec{n} d a$ ⁠.

For the flow in the x-direction, $\vec{u} = u (z, t) \vec{i}, \vec{n} = \vec{i}, d a = d x d z, x \in [0,1], z \in [0,1]$ ⁠, the flow rate is given by $Q_{x} = \iint_{A} u (z, t) d x d z = \int_{0}^{L} \int_{0}^{1} u (z, t) d x d z = L \int_{0}^{1} u (z, t) d z$ ⁠. 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

\begin{matrix} Q_{x} (t) & = & \int_{0}^{1} u (z, t) d z = \sum_{n = 0}^{\infty} \frac{4 \sqrt{2} t^{α}}{{(2 n + 1)}^{2} π^{2} K^{2}} \\ \times \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} t^{α k} {(b_{n}^{2} + 4)}^{k / 2} cos (k α_{n})}{Γ (1 + α (k + 1))} . \end{matrix}

(30)

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

\begin{matrix} Q_{y} (t) & = & \int_{0}^{1} v (z, t) d z = \sum_{n = 0}^{\infty} \frac{4 \sqrt{2} t^{α}}{{(2 n + 1)}^{2} π^{2} K^{2}} \\ \times \sum_{k = 0}^{\infty} \frac{{(- 1)}^{k} t^{α k} {(b_{n}^{2} + 4)}^{k / 2} \sin (k α_{n})}{Γ (1 + α (k + 1))}, \end{matrix}

(31)

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

\begin{matrix} Q_{x} (t) & = & \frac{\sqrt{2}}{64 d^{3} [cos h^{2} (d) - si n^{2} (d)]} \{sinh (d) cosh (d) (1 - cos (2 d))) \\ (- sin (d) cos (d) (1 + cosh (2 d)) + sin h^{2} (d) sin (2 d)) \\ (+ co s^{2} (d) sinh (2 d)\} + \frac{4 \sqrt{2}}{K^{2}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{2} π^{2} (b_{n}^{2} + 4)} 2 sin (2 t) \\ - b_{n} cos (2 t)]) e^{- b_{n} t}, \end{matrix}

(32)

\begin{matrix} Q_{y} (t) & = & \frac{\sqrt{2}}{64 d^{3} [cos h^{2} (d) - si n^{2} (d)]} \{sinh (d) cosh (d) (1 + cos (2 d))) \\ - sin (d) c o s (d) (1 - cosh (2 d)) + (sin h^{2} (d) (1 + \sinh (2 d))) \\ (- \cosh^{2} (d)\} + \frac{4 \sqrt{2}}{K^{2}} \sum_{n = 0}^{\infty} \frac{1}{{(2 n + 1)}^{2} π^{2} (b_{n}^{2} + 4)} [) b_{n} \sin (2 t) \\ + 2 cos (2 t)]) e^{- b_{n} t} . \end{matrix}

(33)

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,

K^{2} (s^{α} + 2 i) \bar{φ} (z, s) = \frac{1}{s} + \frac{\partial^{2} \bar{φ} (z, s)}{\partial z^{2}},

(34)

\bar{φ} (0, s) = \bar{φ} (1, s) = 0,

(35)

whose solution is

\bar{φ} (z, s) = \frac{1}{K^{2} s (s^{α} + 2 i)} \{1 - \frac{cosh [\frac{2 z - 1}{2} \sqrt{K^{2} (s^{α} + 2 i)}]}{\cosh [\frac{1}{2} \sqrt{K^{2} (s^{α} + 2 i)}]}\} .

(36)

It is possible to obtain the inverse Laplace transform of function $\bar{φ} (z, s)$ given by Eq. (36) but the closed form of the corresponding complex velocity field $φ (z, t) = L^{- 1} \{\bar{φ} (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 $\bar{φ} (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²⁸). Due to this reason, in the present paper, we use Stehfest’s algorithm²⁹ for the inverse Laplace transform.

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

φ (z, t) ≅ \frac{ln (2)}{t} \sum_{j = 1}^{2 k} d_{j} \bar{φ} (z, j \frac{\ln (2)}{t}),

(37)

where k is a positive integer,

d_{j} = \sum_{m = [\frac{j + 1}{2}]}^{min (j, k)} \frac{{(- 1)}^{j + k} m^{k} (2 m)!}{(k - m)! m! (m - 1)! (j - m)! (2 m - j)!},

(38)

and $[\frac{j + 1}{2}]$ is the integer part of the number (j + 1)/2.

A comparison of numerical results obtained with Eqs. (15) and (37) will provide information about the accuracy of the numerical calculations.

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.,³⁰ on the computation of the Mittag-Leffler functions. In agreement with the proposed algorithm²⁸ [Eq. (19)], we shall use for the numerical calculations of the values of the Mittag-Leffler function E_α,α+1(z) the formula

E_{α, α + 1} (z) = \frac{1}{z} - \frac{sin (α π)}{α π} \int_{0}^{\infty} \frac{e^{- r^{1 / α}}}{r^{2} - 2 r z cos (α π) + z^{2}} d r, α \in (0,1) .

(39)

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² = 5, α = 0.5a, and p = 10, and a good agreement between the numerical values is found.

TABLE I.

Numerical values of u(z, t) obtained with the formulae (15), (39), and (37).

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² = 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.

FIG. 2.

$FIG. 2. The velocity sketches u(z, t) for t = 0.5 (a), t = 1.5 (b), and t = 3 (c) and various values of the fractional parameter α.$

View large Download slide

The velocity sketches u(z, t) for t = 0.5 (a), t = 1.5 (b), and t = 3 (c) and various values of the fractional parameter α.

FIG. 3.

$FIG. 3. The velocity sketches v(z, t) for t = 0.5 (a), t = 1.5 (b), and t = 3 (c) and various values of the fractional parameter α.$

View large Download slide

The velocity sketches v(z, t) for t = 0.5 (a), t = 1.5 (b), and t = 3 (c) and various values of the fractional parameter α.

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.

FIG. 4.

View large Download slide

The flow rate sketches Q_x(t) for small values of the rotation parameter for α = 0.6 (a) and α = 1.0 (b).

FIG. 5.

View large Download slide

The flow rate profiles Q_y(t) for small values of the rotation parameter for α = 0.6 (a) and α = 1.0 (b).

FIG. 6.

View large Download slide

The flow rate profiles Q_x(t) for large values of the rotation parameter for α = 0.6 (a) and α = 1.0 (b).

FIG. 7.

View large Download slide

The flow rate profiles Q_y(t) for higher values of the rotation parameter for α = 0.6 (a) and α = 1.0 (b).

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.

FIG. 8.

View large Download slide

Comparison of the steady-state solutions given by Eqs. (23) and (28) and Eqs. (26) and (29), respectively, for t = 0.5 (a) and t = 1.0 (b).

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.

FIG. 9.

View large Download slide

Profiles of the velocity u(y, t) for the rotating channel (K = 2) and fixed channel (K = 0) for t = 0.5 (a) and t = 1.0 (b).

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.

FIG. 10.

$FIG. 10. Profiles of the flow rate Qx(t) for a fixed channel (K = 0) and different values of the fractional parameter.$

View large Download slide

Profiles of the flow rate Q_x(t) for a fixed channel (K = 0) and different values of the fractional parameter.

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,

\begin{matrix} \bar{φ} (z, s) & = & \frac{1}{K^{2} s (s^{α} + 2 i)} \\ \times [1 - \frac{sinh (z K \sqrt{s^{α} + 2 i}) + sinh ((1 - z) K \sqrt{s^{α} + 2 i})}{sinh (K \sqrt{s^{α} + 2 i})}] . \end{matrix}

(40)

The comparison between profiles of velocity u(y, t) given by Eqs. (15) and (40) is presented in Fig. 11. It is observed from Fig. 11 that the two solutions are identical, so the accuracy of the series solution is verified.

FIG. 11.

View large Download slide

Comparison between profiles of velocity u(u, t) given by Eq. (15) and Eq. (40), respectively, for t = 0.5 (a) and t = 1.0 (b).

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²⁸ 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 $φ$ (z, t)

Equation (36) can be written as

\bar{φ} (z, s) = {\bar{φ}}_{1} (s) - \frac{1}{s} {\bar{φ}}_{2} (z, s),

(A1)

where

{\bar{φ}}_{1} (s) = \frac{1}{K^{2} s (s^{α} + 2 i)} = \frac{1}{K^{2}} \frac{s^{- 1}}{s^{α} + 2 i} = \frac{1}{K^{2}} \frac{s^{α - (α + 1)}}{s^{α} + 2 i},

(A2)

and

{\bar{φ}}_{2} (z, s) = \frac{1}{w (s)} \frac{cosh [\frac{2 z - 1}{2} \sqrt{w (s)}]}{cosh [\frac{1}{2} \sqrt{w (s)}]}, w (s) = K^{2} (s^{α} + 2 i) .

(A3)

The inverse Laplace transform of function ${\bar{φ}}_{1} (s)$ is

φ_{1} (t) = \frac{1}{K^{2}} t^{α} E_{α, α + 1} (- 2 i t^{α}),

(A4)

where $E_{α, β} (z) = \sum_{n = 0}^{\infty} \frac{z^{n}}{Γ (α n + β)}, L^{- 1} \{\frac{s^{α - β}}{s^{α} - c}\} = t^{β - 1} E_{α, β} (c t^{α})$ denotes the two-parameter Mittag-Leffler function.

To determine the inverse Laplace transform of function ${\bar{φ}}_{2} (z, s)$ ⁠, we consider the auxiliary function

\begin{matrix} \bar{ϕ} (z, s) & = & \frac{cosh (\frac{2 z - 1}{2} \sqrt{s})}{s cosh (\frac{1}{2} \sqrt{s})} = \frac{1}{s} [exp ((z - 1) \sqrt{s})) \\ (- exp (- (z + 1) \sqrt{s})] \frac{1}{1 + exp (- \sqrt{s})} \\ = & \sum_{k = 0}^{\infty} {(- 1)}^{k} [\frac{exp (- (k + 1 - z) \sqrt{s})}{s} + \frac{exp ((k + 1 + z) \sqrt{s})}{s}], \end{matrix}

(A5)

whose inverse Laplace transform is

ϕ (z, t) = \sum_{k = 0}^{\infty} {(- 1)}^{k} [e r f c (\frac{k + 1 - z}{2 \sqrt{t}}) + e r f c (\frac{k + 1 + z}{2 \sqrt{t}})] .

(A6)

Note that ${\bar{φ}}_{2} (z, s) = \bar{ϕ} (z, w (s))$ ⁠, therefore, the inverse Laplace transform of ${\bar{φ}}_{2} (z, s)$ is given by

L^{- 1} \{{\bar{φ}}_{2} (z, s)\} = φ_{2} (z, t) = \int_{0}^{\infty} ϕ (z, x) θ (x, t) d x,

(A7)

where

\begin{matrix} θ (x, t) & = & L^{- 1} \{e^{- x w (s)}\} = L^{- 1} \{exp (- 2 i K^{2} x) exp (- K^{2} x s^{α})\} \\ = & exp (- 2 i K^{2} x) t^{- 1} Θ (0, - α, - K^{2} x s^{- α}), 0 < α < 1 . \end{matrix}

(A8)

In the above relation, Θ(β, −α, −K²xs^−α) is the Wright function given by

L^{- 1} \{s^{- β} exp (- a s^{α})\} = t^{β - 1} Θ (β, - α, - a t^{- α}), 0 < α < 1 .

(A9)

Finally, the inverse Laplace transform of function $\bar{φ} (z, s)$ is given by

φ (z, t) = φ_{1} (t) - \int_{0}^{t} φ_{2} (z, τ) d τ .

(A10)

REFERENCES

1.

C. Y.

Wang

, “

Starting flow in a rotating parallel plate channel

,”

Acta Mech. Sin.

28

(

5

),

1271

–

1276

(

2012

).

https://doi.org/10.1007/s10409-012-0096-5

Google Scholar

Crossref

2.

M. A.

Imran

,

D.

Vieru

, and

I. A.

Mirza

, “

The Influence of Ekman number on flows over an oscillating isothermal vertical plate in a rotating frame

,”

J. Appl. Fluid Mech.

8

(

4

),

781

–

791

(

2015

).

https://doi.org/10.18869/acadpub.jafm.67.223.22941

Google Scholar

Crossref

3.

B. K.

Das

,

M.

Guria

, and

R. N.

Jana

, “

Unsteady Couette flow in a rotating system

,”

Meccanica

43

,

517

–

521

(

2008

).

https://doi.org/10.1007/s11012-008-9130-x

Google Scholar

Crossref

4.

T.

Hayat

,

C.

Fetecau

, and

M.

Sajid

, “

On MHD transient flow of a Maxwell fluid in a porous medium and rotating frame

,”

Phys. Lett. A

372

,

1639

–

1644

(

2008

).

https://doi.org/10.1016/j.physleta.2007.10.036

Google Scholar

Crossref

5.

C. C.

Chang

and

C. Y.

Wang

, “

Rotating electro-osmotic flow over a plate or between two plates

,”

Phys. Rev. E

84

,

056320

(

2011

).

https://doi.org/10.1103/physreve.84.056320

Google Scholar

Crossref

6.

M. V.

Krishna

,

E.

Neeraja

, and

G. S.

Babu

, “

Hall effects on MHD convection flow through a medium in a rotating parallel plate channel

,”

Asian J. Curr. Eng. Math.

2

(

3

),

181

–

189

(

2013

).

Google Scholar

7.

A. M.

Alghamdi

,

S.

Gala

, and

M. A.

Ragusa

, “

A regularity criterion of smooth solution for the 3D viscous Hall-MHD equations

,”

AIMS Math.

3

(

4

),

565

–

574

(

2018

).

https://doi.org/10.3934/Math.2018.4.565

Google Scholar

Crossref

8.

A. M.

Alghamdi

,

S.

Gala

, and

M. A.

Ragusa

, “

New regularity criteria for the 3D Hall-MHD equations

,”

Ann. Polonici Math.

121

(

1

),

7

–

20

(

2018

).

https://doi.org/10.4064/ap170823-23-1

Google Scholar

Crossref

9.

A. M.

Alghamdi

,

S.

Gala

, and

M. A.

Ragusa

, “

On the blow-up criterion for incompressible stokes-MHD equations

,”

Res. Math.

73

(

3

),

110

(

2018

).

https://doi.org/10.1007/s00025-018-0874-x

Google Scholar

Crossref

10.

I. B.

Omrane

,

S.

Gala

,

J.-M.

Kim

, and

M. A.

Ragusa

, “

Logarithmically improved blow-up criterion for smooth solutions to the Leray- α-magnetohydrodynamic equations

,”

Arch. Math.

55

(

1

),

55

–

68

(

2019

).

https://doi.org/10.5817/AM2019-1-55

Google Scholar

Crossref

11.

G. S.

Seth

,

S. M.

Hussain

, and

J. K.

Singh

, “

MHD Couette flow of class-II in a rotating system

,”

J. Appl. Math. Bioinf.

1

(

1

),

31

–

54

(

2011

).

Google Scholar

12.

A.

Razani

, “

Chapman–Jouguet traveling wave for a two-steps reaction scheme

,”

Ital. J. Pure Appl. Math.

39

,

544

–

553

(

2018

).

Google Scholar

13.

E.

Simon

and

N.

Barbara

, “

Schlichtin André Gradient flow formulation and longtime behaviour of a constrained Fokker–Planck equation

,”

Nonlinear Anal.

158

,

142

–

167

(

2017

).

https://doi.org/10.1016/j.na.2017.04.009

Google Scholar

Crossref

14.

S.

Marco

and

W.

Tatsuya

, “

Uniqueness of limit flow for a class of quasi-linear parabolic equations

,”

Adv. Nonlinear Anal.

6

(

2

),

243

–

276

(

2017

).

https://doi.org/10.1515/anona-2016-0134

Google Scholar

Crossref

15.

J.

Simsen

,

M. S.

Simsen

, and

A.

Zimmermann

, “

Aleksandra study of ODE limit problems for reaction-diffusion equations

,”

Opusc. Math.

38

(

1

),

117

–

131

(

2018

).

https://doi.org/10.7494/opmath.2018.38.1.117

Google Scholar

Crossref

16.

D.

Baleanu

,

Fractional Calculus Models and Numerical Methods

(

World Scientific Publishing Company

,

2012

).

Google Scholar

Crossref

17.

I.

Petras

,

Fractional-Order Nonlinear Systems: Modeling Analysis and Simulation

(

Springer and Beijing: HEP

,

London

,

2011

).

Google Scholar

Crossref

18.

A. A.

Kilbas

,

H. M.

Srivastava

, and

J. J.

Trujillo

,

Theory and Applications of Fractional Differential Equations

(

Elsevier

,

Amsterdam, Netherlands

,

2006

).

Google Scholar

19.

L.

Zheng

,

Y.

Liu

, and

X.

Zhang

, “

Slip effects on MHD flow of a generalized Oldroyd-B fluid with fractional derivative

,”

Nonlinear Anal.: Real World Appl.

13

,

513

–

523

(

2012

).

https://doi.org/10.1016/j.nonrwa.2011.02.016

Google Scholar

Crossref

20.

D.

Vieru

,

C.

Fetecau

, and

F.

Corina

, “

Time-fractional free convection flow near a vertical plate with Newtonian heating and mass diffusion

,”

Therm. Sci.

19

(

1

),

85

–

98

(

2015

).

https://doi.org/10.2298/tsci15s1s85v

Google Scholar

Crossref

21.

X. J.

Yang

,

D.

Baleanu

, and

H. M.

Srivastava

,

Local Fractional Integral Transforms and Their Applications

(

Elsevier

,

2015

).

Google Scholar

22.

Y. Z.

Povstenko

, “

Fundamental solutions to time-fractional advection diffusion equation in a case of two space variables

,”

Math. Probl. Eng.

2014

,

705364

.

https://doi.org/10.1155/2014/705364

Crossref

23.

S. K.

Choi

,

B.

Kang

, and

N.

Koo

, “

Stability for Caputo fractional differential systems

,”

Abstr. Appl. Anal.

2014

,

631419

.

https://doi.org/10.1155/2014/631419

Crossref

24.

F.

Zhang

and

Ch.

Li

, “

Stability analysis of fractional differential systems with order lying in (1, 2)

,”

Adv. Differ. Equations

2011

,

213485

https://doi.org/10.1155/2011/213485

Crossref

25.

F.

Mainardi

and

R.

Gorenflo

, “

On Mittag-Leffler type functions in fractional evolution processes

,”

J. Comput. Appl. Math.

118

,

283

–

299

(

2000

).

https://doi.org/10.1016/s0377-0427(00)00294-6

Google Scholar

Crossref

26.

O.

Calin

,

D. C.

Chang

,

K.

Furutani

, and

C.

Iwasaki

, “

The eigenfunction expansion method

,”

Heat Kernels for Eliptic and Sub-Eliptic Operators, Applied and Numerical Harmonic Analysis

(

Birkhauser Boston

,

2010

).

Google Scholar

Crossref

27.

P.

Biswas

and

S.

Singh

, “

Orthogonal eigenfunction expansion method for one-dimensional dual-phase lag heat conduction problem with time-dependent boundary conditions

,”

J. Heat Transfer

140

(

3

),

034501

(

2018

).

https://doi.org/10.1115/1.4037874

Google Scholar

Crossref

28.

B.

Dingfelder

and

J. A. C.

Weideman

, “

An improved Talbot method for numerical Laplace transform inversion

,”

Numer. Algorithms

68

,

167

–

183

(

2015

).

https://doi.org/10.1007/s11075-014-9895-z

Google Scholar

Crossref

29.

H.

Stehfest

, “

Algorithm 368: Numerical inversion of Laplace transform

,”

Commun. ACM

13

(

1

),

47

–

49

(

1970

).

https://doi.org/10.1145/361953.361969

Google Scholar

Crossref

30.

R.

Gorenflo

,

J.

Loutchko

, and

Y.

Luchko

, “

Computation of the Mittag-Leffler function E_α,β(z) and its derivative

,”

Fract. Calc. Appl. Anal.

5

(

4

),

491

–

518

(

2002

).

Google Scholar

2020

Author(s)

Numerical simulations of unsteady flows in a rotating channel using a novel eigenfunction expansion method

I. INTRODUCTION

II. PROBLEM STATEMENT

III. SOLUTION OF THE PROBLEM

A. Eigenfunction expansion of the velocity field

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

C. Flow rates

D. Another form of the complex velocity. Numerical approximation

IV. NUMERICAL RESULTS AND DISCUSSION

V. CONCLUDING REMARKS

DATA AVAILABILITY

ACKNOWLEDGMENTS

APPENDIX: THE SOLUTION OF $φ$ (z, t)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Numerical simulations of unsteady flows in a rotating channel using a novel eigenfunction expansion method

I. INTRODUCTION

II. PROBLEM STATEMENT

III. SOLUTION OF THE PROBLEM

A. Eigenfunction expansion of the velocity field

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

C. Flow rates

D. Another form of the complex velocity. Numerical approximation

IV. NUMERICAL RESULTS AND DISCUSSION

V. CONCLUDING REMARKS

DATA AVAILABILITY

ACKNOWLEDGMENTS

APPENDIX: THE SOLUTION OF φ(z, t)

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX: THE SOLUTION OF $φ$ (z, t)