Axisymmetric Slow Motion of a Porous Spherical Particle in a Viscous Fluid Using Time Fractional Navier–Stokes Equation

Abstract

In this paper, we present the unsteady translational motion of a porous spherical particle in an incompressible viscous fluid. In this case, the modified Navier–Stokes equation with fractional order time derivative is used for conservation of momentum external to the particle whereas modified Brinkman equation with fractional order time derivative is used internal to the particle to govern the fluid flow. Stress jump condition for the tangential stress along with continuity of normal stress and continuity of velocity vectors is used at the porous–liquid interface. The integral Laplace transform technique is employed to solve the governing equations in fluid and porous regions. Numerical inversion code in MATLAB is used to obtain the solution of the problem in the physical domain. Drag force experienced by the particle is obtained. The numerical results have been discussed with the aid of graphs for some specific flows, namely damping oscillation, sine oscillation and sudden motion. Our result shows a significant contribution of the jump coefficient and the fractional order parameter to the drag force.

1. Introduction

Fractional differential equations are a type of differential equation where derivatives are not the traditional derivatives but are of fractional order. Fractional differential equations have several applications in various branches of science and engineering. For example, in a porous media, when there is a fluid flow through it, there is a change in both solid and fluid properties of the porous media due to chemical reactions, mineral precipitation, etc., and, this results in a change in permeability of the porous media and viscosity of the fluid flowing through it over time. The phenomenon that solid and fluid properties change over time is represented by the term ‘memory’. To quantify the effect of history, ‘memory’ is incorporated in the mathematical model. Two types of memory, time memory and space memory, are found in literature. Space memory considers the previous space that the fluids have passed through [1]. One way to include history in any mathematical model is to use fractional order derivative. History of any parameter can be taken into consideration using fractional order derivatives of that parameter. To consider time memory, fractional order derivatives in time are used, and to consider space memory, fractional order derivatives in space are used. Therefore, fractional differential equations are used to model mass transport from fractures in porous media to a porous rock matrix [2]. Fractional differential equations have been used in other different areas such as fractional telegraph equation. The telegraph equation, also known as a damped wave equation is classified as a hyperbolic partial differential equation, which governs physically the voltage and current in an electrical transmission line with distance and time. The fractional telegraph equation has been solved using a combination of homotopy analysis and Laplace transform methods [3]. Similar to the fractional telegraph equation, the equation governing the fluid flow popularly known as the Navier–Stokes equation has been solved for the fractional order time derivative [4]. Kumar et al. [4] developed a new analytical and approximate method by coupling Adomian decomposition method and Laplace transform which they termed as modified Laplace decomposition method to obtain the solution of the fractional Navier–Stokes equation. In their study they replaced the unsteady time derivative term by a fractional order time derivative in the unsteady Navier–Stokes equation. Existence and uniqueness of solutions of the Navier–Stokes equations with fractional time derivatives was studied by Zhou and Peng [5] recently. Wang and Liu [6] developed a modified reduced differential transform method and a new iterative Elzaki transform method to obtain the analytical solutions of the time fractional Navier–Stokes equation while studying the problem of one-dimensional unsteady flow of a viscous fluid in a tube.

The classical Navier–Stokes equation under the assumption of vanishing inertia has been used by many authors to investigate the unsteady motion of spherical objects in viscous fluids. Stokes [7] obtained the drag force on a rigid sphere that is oscillating along one of its diameters in a quiescent fluid. Basset [8] studied the problem of a solid sphere in a viscous fluid undergoing arbitrary translational motion and obtained the drag force. For detailed discussion on the unsteady motion of solid bodies in creeping flows, one may refer to [9] and the references therein. Michaelides and Feng [10] discussed the motion of a small viscous sphere in an unsteady viscous fluid. They reported the presence of two length scales in the momentum diffusion; furthermore, they obtained the force on the sphere and reported that both the length scales appear in the force expression. Choudhuri and Padmavati [11] analyzed the problem of arbitrary unsteady Stokes flow past a liquid sphere. They assumed the shape of the liquid drop to be spherical as a first approximation and used an iterative procedure to determine the shape of the deformed liquid sphere. Bogoyavlenskij [12] obtained smooth and bounded solutions with infinite kinetic energy for the unsteady Navier–Stokes equations of a magneto-hydrodynamic viscous fluid. Ardekani and Rangel [13] studied the unsteady motion of two solid spherical particles in an unbounded incompressible Newtonian flow. They obtained the unsteady force exerted on each particle and reported a larger Besset force due to motion of two spheres rather than the Besset force for a single sphere. Ashmawy [14] investigated the problem of unsteady rotational motion of a rigid sphere about its diameter in a viscous fluid with slip condition. Ashmawy [15] further studied the rotary oscillation of a composite sphere in a viscous liquid bounded by a concentric spherical shell using Navier–Stokes equations and Brinkman model. It was reported that the couple which is responsible for the torque acting on the porous surface decreases with increase in permeability.

Jaber and Ahmad [16] developed analytical solution of the time fractional Navier–Stokes equation using a residual power series method. This power series method is different from the classical power series method in the sense that in the residual power series method one need not compare the coefficients of the corresponding term. They solved the two dimensional problem with variable pressure and obtained the solutions in terms of rapidly convergent series. Their results were found to be in excellent agreement with the existing results. Ashmawy [17] investigated the translational motion of a slip sphere with time-dependent velocity in an incompressible viscous fluid. The unsteady Navier–Stokes equation was used where in the unsteady term the ordinary time derivative was replaced with fractional order time derivative. He obtained an exact formula for the drag force exerted by the fluid on the spherical object and applied this formula to some flows, namely damping oscillation, sine oscillation and sudden motion. It was reported that the order of the fractional derivative contributes considerably to the drag force. An increase in this parameter resulted in an increase in the drag force. Recently, Xu et al. [18] studied the problem of pressure driven flow between two parallel stationary plates. They considered the space fractional Navier–Stokes equation which contains a Riesz space fractional derivative. Their results showed a strong influence of fractional order parameter on the flow velocity and hence, on the flow rate.

There are theoretical and experimental studies available in literature which show that modeling via integer-order derivatives does not provide better prediction in comparison to modeling via fractional order derivatives in the real-world problems. The modeling of physical phenomena via fractional-order derivatives is significant for the control theory, pharmacokinetics, electrical engineering, anomalous diffusion, fluids, electromagnetism, heat transfer (see [19] and references therein). Zhou et al. [20] used the energy method to study the Navier–Stokes equations with the time fractional derivative. Such equations can be used to simulate anomalous diffusion in fractal media. Kang et al. [21] developed a sophisticated diffusion model for gas transport in heterogeneous coal matrix. In view of the disadvantages of existing Fickian diffusion-based models, they proposed an anomalous subdiffusion by introducing the fractional time and space derivatives to classical Fickian flux equation. They concluded that the presence of the fractional parameter describes the gas transport in heterogeneous coal matrix accurately. Muzaffar et al. [22] investigated the helical effects for the fractionalized viscoelastic fluid in helically moved cylinder using fractional derivative in which study of Newtonian and non-Newtonian fluids is presented for rotational and oscillating flows of circular pipe. Zafar and Fetecau [23] investigated the flow of viscous fluid using fractional derivative, here they used integral transforms (Laplace and Fourier sine transforms). Nadeem et al. [24] observed the effects of magnetic field for the flow of a second-grade fluid via fractional derivative in presence of radiative heat transfer. Atangana and Baleanu [25] explored the effects of groundwater flow within confined aquifer by implementing the fractional derivative.

Inspired and motivated from the above discussed studies, the aim of the present study was to investigate the unsteady translational motion of a porous spherical particle in an incompressible viscous fluid which is otherwise at rest. The translational motion of the porous sphere perturbs the fluid flow and keeps it alive. Flow fields in clear fluid and porous regions are governed by unsteady Navier–Stokes and unsteady Brinkman equations with the unsteady term being replaced by fractional time derivative. Hereafter, these equations will be called as modified Navier–Stokes and modified Brinkman equations, respectively. These equations will be used together with mass conservation in both the regions. This work differs from the work of Ashmawy [17] where he dealt with the unsteady translational motion of a solid sphere using fractional time derivative in Navier–Stokes equation. A popularly known stress jump boundary condition [26,27] at the porous–liquid interface is used in place of continuity of stress components. Expression for drag force is obtained. This paper is organized as follows. Section 2 deals with mathematical formulation, solution methodology and implementation of feasible boundary conditions. Section 3 deals with finding the drag force exerted on the porous sphere. Results for different cases of the basic flow are presented in Section 4.

2. Mathematical Formulation and Solution Methodology

Let us consider a porous sphere of radius a and permeability k translating with a time-dependent velocity $U\left(t\right)$ in an incompressible viscous fluid which is otherwise at rest. In addition, we shall assume that the velocity $U\left(t\right)$ takes the form $U\left(t\right)=U_{0}f\left(t\right)$, where the constant $U_0$ has the dimensions of velocity. We have transformed the classical unsteady Stokes and Brinkman equations into modified Stokes and Brinkman equations where the unsteady term has been replaced by a fractional time derivative with the help of the Caputo operator (for more details see [17,19,28,29] and references therein). The analytical solution for velocity and pressure of the time fractional Stokes and time fractional Brinkman equations is examined by Laplace transform technique. The flow fields exterior and interior to the porous sphere under the creeping flow conditions are governed by modified Navier–Stokes equation and modified Brinkman equation together with mass conservation, respectively. Hence, the equations governing the fluid flows are given by

$$\rho {\displaystyle \frac{{\partial }^{\alpha }{\mathbf{v}}^e}{\partial t^{\alpha }}}=-\nabla p^e+\mu {\nabla }^2{\mathbf{v}}^e,\mathrm{and}\nabla ·{\mathbf{v}}^e=0\mathrm{if}r>a$$

(1)

$$\rho {\displaystyle \frac{{\partial }^{\alpha }{\mathbf{v}}^i}{\partial t^{\alpha }}}=-\nabla p^i+\mu {\nabla }^2{\mathbf{v}}^i-{\displaystyle \frac{\mu }{k}}{\mathbf{v}}^i,\mathrm{and}\nabla ·{\mathbf{v}}^i=0\mathrm{if}r<a$$

(2)

where $\rho$ is the density of the fluid, ${\mathbf{v}}^j,j=e,i$ denotes the velocity field. The superscript $j=e$ corresponds to the velocity field exterior to the porous sphere and $j=i$ corresponds to the velocity field inside the porous sphere. Similarly, $p^j,j=e,i$ denotes the pressure at any point in the exterior and interior of the porous sphere, $\mu$ is the dynamic viscosity and k is the permeability. The parameter $\alpha$ is the order of Caputo fractional derivative such that $0\le \alpha \le 1$. The fractional derivative of $f\left(t\right)$ in Caputo sense is defined as

$$\frac{{\partial }^{\alpha }f\left(t\right)}{\partial t^{\alpha }}=\frac{1}{\mathsf{\Gamma }(n-\alpha )}{\int }_0^t\frac{f^{\left(n\right)}\left(x\right)}{{(t-x)}^{\alpha -n+1}}\mathrm{d}x,n-1\le \alpha \le n,n\in \mathbb{N},t>0,$$

where $f^{\left(n\right)}$ represents the nth ordinary derivative of $f\left(t\right)$ with respect to t. One of the great advantages of the Caputo fractional derivative is that it allows traditional initial and boundary conditions to be included in the formulation of the problem. In addition its derivative for a constant is zero. To compute the fractional derivative of a function in the Caputo sense, we must first calculate its derivative. Caputo derivatives are defined only for differentiable functions. Since the velocity and pressure fields in the fluid and porous regions are smooth functions, use of Caputo fractional derivative is justified. We use the following dimensionless variables $\tilde{\mathbf{X}}={\displaystyle \frac{\mathbf{X}}{a}},\tilde{\mathbf{v}}={\displaystyle \frac{\mathbf{v}}{{\mathbf{U}}_{\mathbf{0}}}}$, $\tilde{p}={\displaystyle \frac{p}{(\mu U_0/a)}}$, $\tilde{t}={\displaystyle \frac{t}{t_c}}$ and ${\tilde{\psi }}^e={\displaystyle \frac{{\psi }^e}{U_0{a}^2}}$ to non-dimensionalize the governing equations. Here, $U_0$ is the amplitude of the velocity $U\left(t\right)=U_{0}f\left(t\right)$. As a result, Equations (1) and (2) are transformed into non-dimensional form as

$$T_u{\displaystyle \frac{{\partial }^{\alpha }{\tilde{\mathbf{v}}}^e}{\partial {\tilde{t}}^{\alpha }}}=-\tilde{\nabla }{\tilde{p}}^e+{\tilde{\nabla }}^2{\tilde{\mathbf{v}}}^e,\tilde{\nabla }·{\tilde{\mathbf{v}}}^e=0,$$

(3)

$$T_u{\displaystyle \frac{{\partial }^{\alpha }{\tilde{\mathbf{v}}}^i}{\partial {\tilde{t}}^{\alpha }}}=-\tilde{\nabla }{\tilde{p}}^i+{\tilde{\nabla }}^2{\tilde{\mathbf{v}}}^i-{\lambda }^2{\tilde{\mathbf{v}}}^i,\tilde{\nabla }·{\tilde{\mathbf{v}}}^i=0,$$

(4)

where $T_u={\displaystyle \frac{\rho a^2}{\mu t_c}}$ is the unsteadiness parameter, $\tilde{\nabla }=a\nabla ,{\tilde{\nabla }}^2=a^2{\nabla }^2$ and $t_c$ is a characteristic time scale, ${\lambda }^2={\displaystyle \frac{1}{Da}}$ and $Da={\displaystyle \frac{k}{a^2}}$ is the Darcy number. The porous particle is assumed to be nondeformable. The spherical coordinate system $(r,\theta ,\varphi )$ is established with its origin at the particle center. Using the spherical coordinates $(r,\theta ,\varphi )$, and dropping the symbol $\tilde{}$, the r, $\theta$ and $\varphi$ components in Equations (3) and (4) can be written in dimensionless form as

$$\begin{array}{ccc}{T}_u{\displaystyle \frac{{\partial }^{\alpha }{v}_r^e}{\partial t^{\alpha }}}\hfill & =& -{\displaystyle \frac{\partial p^e}{\partial r}}+{\nabla }^2{v}_r^e,T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\theta }^e}{\partial t^{\alpha }}}=-\frac{1}{r}{\displaystyle \frac{\partial p^e}{\partial \theta }}+{\nabla }^2{v}_{\theta }^e,T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\varphi }^e}{\partial t^{\alpha }}}=-\frac{1}{rsin\theta }{\displaystyle \frac{\partial p^e}{\partial \varphi }}+{\nabla }^2{v}_{\varphi }^e\hfill \\ \nabla ·{\mathbf{v}}^e\hfill & =& 0.\hfill \\ T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_r^i}{\partial t^{\alpha }}}\hfill & =& -{\displaystyle \frac{\partial p^i}{\partial r}}+{\nabla }^2{v}_r^i-{\lambda }^2{v}_r^i,T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\theta }^i}{\partial t^{\alpha }}}=-\frac{1}{r}{\displaystyle \frac{\partial p^i}{\partial \theta }}+{\nabla }^2{v}_{\theta }^i-{\lambda }^2{v}_{\theta }^i\hfill \end{array}$$

(5)

$$\begin{array}{ccc}{T}_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\varphi }^e}{\partial t^{\alpha }}}\hfill & =& -\frac{1}{rsin\theta }{\displaystyle \frac{\partial p^i}{\partial \varphi }}+{\nabla }^2{v}_{\varphi }^i-{\lambda }^2{v}_{\varphi }^i,\nabla ·{\mathbf{v}}^i=0.\hfill \end{array}$$

(6)

Since the flow of the fluid is axially symmetric, all the physical quantities are independent of $\varphi$. Hence, we assume the velocity vectors ${\mathbf{v}}^e$ and ${\mathbf{v}}^i$ of the form

$${\mathbf{v}}^j=v_r^j(r,\theta ){\mathbf{e}}_r+v_{\theta }^j(r,\theta ){\mathbf{e}}_{\theta },j=e,i,$$

(7)

where $({\mathbf{e}}_r,{\mathbf{e}}_{\theta },{\mathbf{e}}_{\varphi })$ are unit base vectors in spherical polar coordinates. Introducing the stream functions ${\psi }^j(r,\theta ),j=e,i$, through

$$v_r^j=-{\displaystyle \frac{1}{r^{2}sin\theta }}{\displaystyle \frac{\partial {\psi }^j}{\partial \theta }}\mathrm{and}{v}_{\theta }^j={\displaystyle \frac{1}{rsin\theta }}{\displaystyle \frac{\partial {\psi }^j}{\partial r}},j=e,i$$

(8)

Using Equation (8), Equations (5) and (6) become

$$\begin{array}{ccc}\hfill T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_r^e}{\partial t^{\alpha }}}& =& -{\displaystyle \frac{\partial p^e}{\partial r}}-{\displaystyle \frac{1}{r^{2}sin\theta }}{\displaystyle \frac{\partial }{\partial \theta }}\left(E^2{\psi }^e\right)\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\theta }^e}{\partial t^{\alpha }}}& =& -{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial p^e}{\partial \theta }}+{\displaystyle \frac{1}{rsin\theta }}{\displaystyle \frac{\partial }{\partial r}}\left(E^2{\psi }^e\right)\hfill \end{array}$$

(10)

$$\frac{1}{r^2}\frac{\partial \left(r^2{v}_r^e\right)}{\partial r}+\frac{1}{rsin\theta }\frac{\partial (v_{\theta }^{e}sin\theta )}{\partial \theta }=0$$

(11)

$$\begin{array}{ccc}\hfill T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_r^i}{\partial t^{\alpha }}}& =& -{\displaystyle \frac{\partial p^i}{\partial r}}-{\displaystyle \frac{1}{r^{2}sin\theta }}{\displaystyle \frac{\partial }{\partial \theta }}\left(E^2{\psi }^i\right)+{\displaystyle \frac{{\lambda }^2}{r^{2}sin\theta }}{\displaystyle \frac{\partial {\psi }^i}{\partial \theta }}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill T_u{\displaystyle \frac{{\partial }^{\alpha }{v}_{\theta }^i}{\partial t^{\alpha }}}& =& -{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial p^i}{\partial \theta }}+{\displaystyle \frac{1}{rsin\theta }}{\displaystyle \frac{\partial }{\partial r}}\left(E^2{\psi }^i\right)-{\displaystyle \frac{{\lambda }^2}{rsin\theta }}{\displaystyle \frac{\partial {\psi }^i}{\partial r}}\hfill \end{array}$$

(13)

$$\frac{1}{r^2}\frac{\partial \left(r^2{v}_r^i\right)}{\partial r}+\frac{1}{rsin\theta }\frac{\partial (v_{\theta }^{i}sin\theta )}{\partial \theta }=0$$

(14)

where $E^2$ is a linear operator defined as follows

$$E^2={\displaystyle \frac{{\partial }^2}{\partial r^2}}-{\displaystyle \frac{cot\theta }{r^2}}{\displaystyle \frac{\partial }{\partial \theta }}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{{\partial }^2}{\partial {\theta }^2}}$$

We solve the Equations (9) and (10), and the Equations (12) and (13) using the integral Laplace transform. The Laplace transform of a function $h(r,\theta ,t))$ is defined by

$$\overline{h}(r,\theta ,s)={\int }_0^{\infty }h(r,\theta ,t)e^{-st}dt$$

(15)

$$L\left\{{\displaystyle \frac{{\partial }^{\alpha }h(r,\theta ,t)}{\partial t^{\alpha }}}\right\}=s^{\alpha }\overline{h}(r,\theta ,s)-s^{\alpha -1}h(r,\theta ,0)$$

(16)

Taking the Laplace transform of the Equations (9) and (10), and the Equations (12) and (13), we arrive at the following equations

$$\begin{array}{ccc}\hfill T_u{s}^{\alpha }{\overline{v}}_r^e& =& -{\displaystyle \frac{\partial {\overline{p}}^e}{\partial r}}-{\displaystyle \frac{1}{r^{2}sin\theta }}{\displaystyle \frac{\partial }{\partial \theta }}{E}^2{\overline{\psi }}^e\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill T_u{s}^{\alpha }{\overline{v}}_{\theta }^e& =& -{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial {\overline{p}}^e}{\partial \theta }}+{\displaystyle \frac{1}{rsin\theta }}{\displaystyle \frac{\partial }{\partial r}}{E}^2{\overline{\psi }}^e\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill T_u{s}^{\alpha }{\overline{v}}_r^i& =& -{\displaystyle \frac{\partial {\overline{p}}^i}{\partial r}}-{\displaystyle \frac{1}{r^{2}sin\theta }}{\displaystyle \frac{\partial }{\partial \theta }}{E}^2{\overline{\psi }}^i+{\displaystyle \frac{{\lambda }^2}{r^{2}sin\theta }}{\displaystyle \frac{\partial {\overline{\psi }}^i}{\partial \theta }}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill T_u{s}^{\alpha }{\overline{v}}_{\theta }^i& =& -{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial {\overline{p}}^i}{\partial \theta }}+{\displaystyle \frac{1}{rsin\theta }}{\displaystyle \frac{\partial }{\partial r}}{E}^2{\overline{\psi }}^i-{\displaystyle \frac{{\lambda }^2}{rsin\theta }}{\displaystyle \frac{\partial {\overline{\psi }}^i}{\partial r}}\hfill \end{array}$$

(20)

where the following conditions have been used

$$v_r^e(r,\theta ,0)=v_{\theta }^e(r,\theta ,0)=0,and{v}_r^i(r,\theta ,0)=v_{\theta }^i(r,\theta ,0)=0$$

(21)

Eliminating the pressure form (17) and (18), we get the following equation

$$E^2\left(E^2-l^2\right){\overline{\psi }}^e=0,\mathrm{where}{l}^2=s^{\alpha }{T}_u$$

(22)

similarly, after eliminating the pressure from (19) and (20), we obtain

$$E^2\left(E^2-m^2\right){\overline{\psi }}^i=0,\mathrm{where}{m}^2=(l^2+{\lambda }^2).$$

(23)

The general solution of the partial differential equation (PDE), defined in Equation (22) is given as follows,

$${\overline{\psi }}^e=\left\{{\displaystyle \frac{A_1}{r}}+A_2{r}^2+A_3\sqrt{r}{I}_{3/2}\left(lr\right)+A_4\sqrt{r}{K}_{3/2}\left(lr\right)\right\}{sin}^2\theta ,$$

(24)

where $I_n(.)$ and $K_n(.)$ are the modified Bessel functions of the first and second kinds of order n and $A_1,A_2,A_3,A_4$ are the unknown constants to be determined using feasible boundary conditions at the porous -liquid interface. Since the solution of the problem must be bounded as r tends to infinity, we must have $A_2=A_3=0$. Similarly, the general solution of Equation (23) can be written as,

$${\overline{\psi }}^i=\left\{{\displaystyle \frac{B_1}{r}}+B_2{r}^2+B_3\sqrt{r}{I}_{3/2}\left(mr\right)+B_4\sqrt{r}{K}_{3/2}\left(mr\right)\right\}{sin}^2\theta ,$$

(25)

where $B_1,B_2,B_3,B_4$ are the unknown constants to be determined. Here also the solution must be bounded as $r\to 0$, therefore we set $B_1=B_4=0$. Substituting the stream functions given in Equations (24) and (25) into Equation (8), we get the velocity components external and internal to the porous sphere as follows

$$\begin{array}{ccc}\hfill {\overline{v}}_r^e& =& -2\left\{{\displaystyle \frac{A_1}{r^3}}+A_4{r}^{-3/2}{K}_{3/2}\left(lr\right)\right\}cos\theta ,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\overline{v}}_{\theta }^e& =& -\left\{{\displaystyle \frac{A_1}{r^3}}+A_4{r}^{-3/2}\left(K_{3/2}\left(lr\right)+lr{K}_{1/2}\left(lr\right)\right)\right\}sin\theta \hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill {\overline{v}}_r^i& =& -2\left\{B_2+B_3{r}^{-3/2}{I}_{3/2}\left(mr\right)\right\}cos\theta ,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\overline{v}}_{\theta }^i& =& \left\{2B_2-B_3{r}^{-3/2}\left(I_{3/2}\left(mr\right)-mr{I}_{1/2}\left(mr\right)\right)\right\}sin\theta \hfill \end{array}$$

(29)

Substitution of Equations (24) and (25) into Equations (18) and (20), respectively, and subsequent integration of these equations yield the pressure expressions exterior and interior to the porous sphere which are given as follows

$$\begin{array}{ccc}\hfill {\overline{p}}^e& =& p_0-A_1{l}^2{r}^{-2}cos\theta \hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill {\overline{p}}^i& =& p_0+2B_2{m}^{2}rcos\theta ,\hfill \end{array}$$

(31)

where $p_0$ is a constant.

The following subsection presents a short discussion on the use of boundary conditions at the porous–liquid interface.

Boundary Conditions

The applicability of a suitable boundary condition at the porous–liquid interface is a topic of debate. However, researchers have used a set of boundary conditions depending on the nature of porous material and the corresponding governing equations that are employed inside the porous medium. For example, Darcy’s law is assumed to govern the fluid flow through porous media with low permeability and the same in case of high permeability is governed by generalized Darcy’s (Brinkman) equation. Due to the difference in the order of Darcy equation and Stokes equation, the boundary condition in terms of stresses cannot be imposed in case of Stokes–Darcy coupling. Instead of conditions on stresses, pressure continuity is used along with continuity of normal velocity and slip in tangential velocity (please see [30] and references therein). However, the case of Stokes–Brinkman coupling enables one to use stress boundary conditions due to matching order of Stokes and Brinkman equations, i.e., 2. In general, for a Stokes–Brinkman coupled system, continuity of stresses together with continuity of velocity components are used at the porous–liquid interface. One may refer Nield [31] for developments on the boundary conditions at a porous–liquid interface. For the current problem under investigation, we use continuity of velocity components and continuity of normal stress component together with a jump in tangential stress popularly known as stress jump condition (please see [26,27,32] and references therein for more details). Accordingly, we use the following boundary conditions at a porous–liquid interface in order to determine the velocity and pressure fields outside and inside the porous sphere. In addition, we shall assume that the velocity of the boundary takes the form $U\left(t\right)=U_{0}f\left(t\right)$, where $U_0$ is a constant.

We apply the following boundary conditions on $r=1$:

(i)

continuity of the velocity components: ${\overline{v}}_r^e-{\overline{v}}_r^i=U\left(t\right)cos\theta ,{\overline{v}}_r^e-{\overline{v}}_r^i=-U\left(t\right)sin\theta$,

(ii)

continuity of the normal stress component: ${\overline{T}}_{rr}^e={\overline{T}}_{rr}^i,$

(iii)

stress jump boundary condition for the shear stress:

$$\begin{array}{c}\hfill {\overline{T}}_{r\theta }^e-{\overline{T}}_{r\theta }^i=\beta \lambda {\overline{v}}_{\theta }^i,\end{array}$$

(32)

where $\beta$ is the stress jump coefficient. The unknowns $A_1,A_4$ and $B_2,B_3$ are evaluated with the help of the aforementioned boundary conditions. These coefficients are obtained in Laplace domain and they are given as follows:

$$A_1=\frac{M_1}{M_2},$$

(33)

where

$$\begin{array}{ccc}\hfill M_1& =& -\overline{f}\left(s\right)(-24cosh\left(m\right)\beta \lambda m+21sinh\left(m\right)\beta \lambda m^2-\beta \lambda l^2{m}^{4}sinh\left(m\right)+\beta \lambda l^2{m}^{3}cosh\left(m\right)\hfill \\ & & +7\beta \lambda l^2{m}^{2}sinh\left(m\right)-3cosh\left(m\right)l{m}^5+12sinh\left(m\right)l{m}^2+3cosh\left(m\right)\beta \lambda m^3-\hfill \\ & & 3sinh\left(m\right)\beta \lambda m^4+15sinh\left(m\right)l{m}^4+24sinh\left(m\right)\beta \lambda l-12cosh\left(m\right)l{m}^3+\hfill \\ & & 3l^2{m}^{2}sinh\left(m\right)-8\beta \lambda l^{2}mcosh\left(m\right)-3sinh\left(m\right)\beta \lambda l{m}^4+3cosh\left(m\right)\beta \lambda l{m}^3+\hfill \\ & & 21sinh\left(m\right)\beta \lambda l{m}^2-24cosh\left(m\right)\beta \lambda lm-3cosh\left(m\right)m^5+15sinh\left(m\right)m^4-12cosh\left(m\right)m^3\hfill \\ & & +12sinh\left(m\right)m^2+24sinh\left(m\right)\beta \lambda -l^3{m}^{2}sinh\left(m\right)+8\beta \lambda l^{2}sinh\left(m\right)-3l^2{m}^{3}cosh\left(m\right)+\hfill \\ & & l^3{m}^{3}cosh\left(m\right)-l^2{m}^{5}cosh\left(m\right)+6l^2{m}^{4}sinh\left(m\right)+l^3{m}^{4}sinh\left(m\right))\hfill \end{array}$$

$$\begin{array}{ccc}\hfill M_2& =& l^2(-9cosh\left(m\right)\beta \lambda m+9sinh\left(m\right)\beta \lambda m^2-\beta \lambda l^2{m}^{2}sinh\left(m\right)+43sinh\left(m\right)l{m}^2+\hfill \\ & & 2cosh\left(m\right)\beta \lambda m^3-2sinh\left(m\right)\beta \lambda m^4+2sinh\left(m\right)l{m}^4-7sinh\left(m\right)\beta \lambda l-7cosh\left(m\right)l{m}^3\hfill \\ & & -36cosh\left(m\right)m-l^{3}sinh\left(m\right)+l^{3}mcosh\left(m\right)+6l^2{m}^{2}sinh\left(m\right)+3l^{2}sinh\left(m\right)+36sinh\left(m\right)\hfill \\ & & +3\beta \lambda l^{2}mcosh\left(m\right)-5sinh\left(m\right)\beta \lambda l{m}^2+7cosh\left(m\right)\beta \lambda lm-2cosh\left(m\right)m^5+12sinh\left(m\right)m^4\hfill \\ & & -15cosh\left(m\right)m^3+51sinh\left(m\right)m^2+9sinh\left(m\right)\beta \lambda +36sinh\left(m\right)l+l^3{m}^{2}sinh\left(m\right)-\hfill \\ & & 3\beta \lambda l^{2}sinh\left(m\right)-3l^{2}mcosh\left(m\right)-l^2{m}^{3}cosh\left(m\right)-36cosh\left(m\right)lm)\hfill \end{array}$$

$$A_4=\frac{M_3}{M_4},$$

(34)

where

$$\begin{array}{ccc}\hfill M_3& =& -\sqrt{2}{\mathrm{e}}^l{l}^{3/2}\overline{f}\left(s\right)(\beta \lambda l^2{m}^{2}sinh\left(m\right)-3sinh\left(m\right)\beta \lambda m^4-3cosh\left(m\right)m^5+15sinh\left(m\right)m^4+\hfill \\ & & \beta \lambda l^{2}mcosh\left(m\right)+3cosh\left(m\right)\beta \lambda m^3-\beta \lambda l^{2}sinh\left(m\right)+21sinh\left(m\right)\beta \lambda m^2-12cosh\left(m\right)m^3\hfill \\ & & +12sinh\left(m\right)m^2-24cosh\left(m\right)\beta \lambda m+24sinh\left(m\right)\beta \lambda )\hfill \\ \hfill M_4& =& \sqrt{\pi }{M}_2.\hfill \end{array}$$

$$B_2=\frac{M_5}{M_6},$$

(35)

where

$$\begin{array}{ccc}\hfill M_5& =& \overline{f}\left(s\right)l^2(\beta \lambda l^2{m}^{2}sinh\left(m\right)-l^3{m}^{2}sinh\left(m\right)+l^2{m}^{3}cosh\left(m\right)+3sinh\left(m\right)\beta \lambda l{m}^2-\hfill \\ & & 6l^2{m}^{2}sinh\left(m\right)-\beta \lambda l^{2}mcosh\left(m\right)-l^{3}mcosh\left(m\right)+3cosh\left(m\right)l{m}^3+\beta \lambda l^{2}sinh\left(m\right)+\hfill \\ & & 3sinh\left(m\right)\beta \lambda m^2+l^{3}sinh\left(m\right)-39sinh\left(m\right)l{m}^2-3cosh\left(m\right)\beta \lambda lm+3l^{2}mcosh\left(m\right)+\hfill \\ & & 3cosh\left(m\right)m^3+3sinh\left(m\right)\beta \lambda l-3l^{2}sinh\left(m\right)-39sinh\left(m\right)m^2-3cosh\left(m\right)\beta \lambda m+\hfill \\ & & 36cosh\left(m\right)lm+3sinh\left(m\right)\beta \lambda -36sinh\left(m\right)l+36cosh\left(m\right)m-36sinh\left(m\right))\hfill \end{array}$$

$$\begin{array}{ccc}\hfill M_6& =& 2l^2(9cosh\left(m\right)\beta \lambda m-9sinh\left(m\right)\beta \lambda m^2+\beta \lambda l^2{m}^{2}sinh\left(m\right)-43sinh\left(m\right)l{m}^2-\hfill \\ & & 2cosh\left(m\right)\beta \lambda m^3+2sinh\left(m\right)\beta \lambda m^4-2sinh\left(m\right)l{m}^4+7sinh\left(m\right)\beta \lambda l+7cosh\left(m\right)l{m}^3+\hfill \\ & & 36cosh\left(m\right)m+l^{3}sinh\left(m\right)-l^{3}mcosh\left(m\right)-6l^2{m}^{2}sinh\left(m\right)-3l^{2}sinh\left(m\right)-36sinh\left(m\right)\hfill \\ & & -3\beta \lambda l^{2}mcosh\left(m\right)+5sinh\left(m\right)\beta \lambda l{m}^2-7cosh\left(m\right)\beta \lambda lm+2cosh\left(m\right)m^5-12sinh\left(m\right)m^4\hfill \\ & & +15cosh\left(m\right)m^3-51sinh\left(m\right)m^2-9sinh\left(m\right)\beta \lambda -36sinh\left(m\right)l-l^3{m}^{2}sinh\left(m\right)+\hfill \\ & & 3\beta \lambda l^{2}sinh\left(m\right)+3l^{2}mcosh\left(m\right)+l^2{m}^{3}cosh\left(m\right)+36cosh\left(m\right)lm)\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill B_3& =& \frac{-\overline{f}\left(s\right)m^{3/2}\sqrt{\pi }{l}^2\left(\beta \lambda l^2-3l{m}^2+3\beta \lambda l-3m^2+3\beta \lambda \right)}{{\displaystyle \frac{M_6}{\sqrt{2}}}},\hfill \end{array}$$

(36)

where $\overline{f}\left(s\right)$ denotes the Laplace transform of the function $f\left(t\right)$ and $l,m$ are functions of s in Laplace domain. It is worth mentioning here that to obtain any physical quantity such as force and velocities in the physical domain, we need to obtain the inverse Laplace transform of the coefficients. In the following section, we calculate force exerted by the surrounding fluid on the surface of porous particle.

3. Force Acting on the Particle

The drag force (in the z direction) exerted by the external fluid on the porous sphere with the spherical boundary $r=1$ in dimensionless form in the Laplace domain can be determined from

$$\begin{array}{ccc}\hfill \overline{F_z^{\alpha }}& =& {\int }_{\varphi =0}^{2\pi }{\int }_{\theta =0}^{\pi }\left[T_{rr}^e\widehat{{\mathbf{e}}_r}+T_{r\theta }^e\widehat{{\mathbf{e}}_{\theta }}+T_{r\varphi }^e\widehat{{\mathbf{e}}_{\varphi }}\right]r^{2}sin\theta d\theta d\varphi {|}_{r=1}\hfill \\ & =& \frac{4\pi }{3}{l}^2\left[A_1-2A_4{K}_{3/2}\left(l\right)\right].\hfill \end{array}$$

(37)

We take an inverse Laplace transform of $A_1$ and $A_4$ to get the force in the time domain. It may be noted here that $Da\to 0$ or $\lambda \to \infty$ corresponds to the case of solid sphere and force obtained in Equation (37) reduces to

$$\overline{F_z^{\alpha }}=-6\pi \left(\overline{f}\left(s\right)+\frac{1}{9}{l}^2\overline{f}\left(s\right)+l\overline{f}\left(s\right)\right).$$

(38)

It may be noted that the expression obtained in Equation (38) in dimensionless form agrees with Equation (18) as a limiting case of $\gamma \to \infty$ in [17]. In this case, the force can be obtained in the time domain by following the procedure mentioned in [17]. Moreover, the classical formula of the drag force in the case of no-slip given by Basset [8] can be deduced as a special case by assuming $Da\to 0$ and $\alpha =1$ in Equation (37) and taking the inverse Laplace transform as stated in [17]. However, for the general case, we obtain a cumbersome force expression in Laplace domain and due to inability to find out closed form expression after taking inverse Laplace transform, we resort to numerical inverse Laplace transform. The Gaver–Stehfest algorithm [33] has so far been the most popular technique to compute the inverse Laplace transform in the context of transient electromagnetics. However, the accuracy of the Gaver–Stehfest algorithm, even when using double-precision arithmetic, is relatively low at late times due to round-off errors. In addition, the Gaver–Stehfest algorithm is significantly problem dependent. Therefore, we have turned our attention to two other algorithms for computing inverse Laplace transforms, namely, the Talbot and Euler algorithms (please see [34,35]). These two methods are the most widely used methods for numerical inversion of Laplace transform due to their efficiency and computational cost. Euler and Talbot have same efficiency; however, Euler has less CPU time compared to Talbot. Therefore, we use the Euler scheme to find inverse Laplace transform numerically. Here we present the Euler algorithm in brief.

Numerical Inverse Laplace Transform (Euler’s Method)

Let f be non-negative real valued and $\overline{F}$ be complex valued functions. Then the Laplace transform of f is defined as,

$$\overline{F}\left(s\right)={\int }_0^{\infty }{e}^{-st}f\left(t\right)dt,\mathrm{where}s\in \mathbb{C}.$$

(39)

Under some sufficient conditions on f the transform (39) is well defined. The inverse transform of $\overline{F}$ is defined in the form of Bromwich inversion integral,

$$f\left(t\right)=\underset{y\to \infty }{lim}{\displaystyle \frac{1}{2\pi i}}{\int }_{a-iy}^{a+iy}{e}^{st}\overline{F}\left(s\right)ds,$$

(40)

where $a\in \mathbb{R}$ and $\overline{F}$ has no singular points on the right side of the vertical line $s=a$. After change of variables, Equation (40) is written in the following form [36]

$$f\left(t\right)=\frac{2}{\pi }{e}^{at}{\int }_0^{\infty }Re\left(\overline{F}(a+iu)\right)cos\left(ut\right)du,\mathrm{where}i=\sqrt{-1}.$$

(41)

There exist many [37,38,39] numerical schemes to find inversion of $\overline{F}$, defined in the Equation (39). However, it has been noticed that the Euler’s method [38,40,41] is computationally efficient as compared to other numerical schemes. Therefore we discuss Euler’s method briefly in the following lines for the sake of continuity.

In Euler’s method, the improper integral defined in the Equation (41) is evaluated using the trapezoidal rule as follows

$$f\left(t\right)\approx f_h\left(t\right)={\displaystyle \frac{h}{\pi }}{e}^{at}Re\left(\overline{F}\right)\left(a\right)+{\displaystyle \frac{2h}{\pi }}{e}^{at}\sum _{k=0}^{\infty }Re\left(\overline{F}\right)(a+ikh)cos\left(kht\right).$$

(42)

Equation (42) can be written in the form of alternating series with the step size $h={\displaystyle \frac{\pi }{2t}}$ and $a={\displaystyle \frac{A}{2t}}$ as follows

$$f_h\left(t\right)={\displaystyle \frac{e^{A/2}}{2t}}Re\left(\overline{F}\right)\left({\displaystyle \frac{A}{2t}}\right)+{\displaystyle \frac{e^{A/2}}{t}}\sum _{k=0}^{\infty }{(-1)}^{k}Re\left(\overline{F}\right)\left({\displaystyle \frac{A+2k\pi i}{2t}}\right).$$

(43)

Equation (43) has been evaluated efficiently using the following Euler summation

$$E(m,n,t)=\sum _{k=0}^m{2}^{-m}\left(\genfrac{}{}{0pt}{}{m}{k}\right)s_{n+k}\left(t\right),$$

(44)

which is defined as the weighted average of the last m partial sums by binomial probability distribution with probability $p={\displaystyle \frac{1}{2}}$ and the partial sums are defined as

$$s_n\left(t\right)={\displaystyle \frac{e^{A/2}}{2t}}Re\left(\overline{F}\right)\left({\displaystyle \frac{A}{2t}}\right)+{\displaystyle \frac{e^{A/2}}{t}}\sum _{k=0}^n{(-1)}^{k}Re\left(\overline{F}\right)\left({\displaystyle \frac{A+2k\pi i}{2t}}\right).$$

(45)

Assuming $n=2m$ in Equation (44), an alternative formula to find the Laplace inverse transform numerically is proposed in the unified framework [35]

$$f\left(t\right)\approx f_n\left(t\right)\equiv {\displaystyle \frac{1}{t}}\sum _{k=0}^{2m}{w}_k\overline{F}\left({\displaystyle \frac{{\alpha }_k}{t}}\right),t>0,$$

(46)

where $w_k$ and ${\alpha }_k$ are the weights and nodes, respectively. The weights and nodes are defined as $w_k={10}^{m/3}{\eta }_k$, ${\eta }_k={(-1)}^k{\xi }_k$, ${\alpha }_k={\displaystyle \frac{mln\left(10\right)}{3}}+\pi ik$, ${\xi }_0={\displaystyle \frac{1}{2}},{\xi }_k=1,1\le k\le m,{\xi }_{2m}={\displaystyle \frac{1}{2^m}}$, and

$${\xi }_{2m-k}={\xi }_{2m-k+1}+2^{-m}\left(\genfrac{}{}{0pt}{}{m}{k}\right),0<k<m.$$

Here $w_k$ and ${\alpha }_k$ are called the exterior and interior scaling factors, they depend only on m and not on the transform $\overline{F}$ or the time argument t. However, the accuracy of numerical scheme, defined in Equation (46), depends on $\overline{F}$ and t.

4. Special Fluid Flows

Here, we present the non-dimensional drag force ${\tilde{F}}_z^{\alpha }=-{\overline{F}}_z^{\alpha }/6\pi \mu U_{0}a$ exerted on the porous sphere for three different types of fluid flow:

Case 1: Damping oscillation

As already mentioned, he velocity of the boundary of the porous sphere is of the form $U\left(t\right)=U_{0}f\left(t\right)$. In this case, we assume that the porous sphere starts to move with the velocity $U\left(t\right)$, but with $f\left(t\right)=e^{-\omega t}sin\omega t$, where $\omega$ is the frequency of oscillations. It may be noted from Equation (37) that the drag force is obtained in the Laplace domain and it depends on the coefficients $A_1$ and $A_4$. It may also be noted that these coefficients are obtained in Equations (33) and (34) which contain Laplace transform of $f\left(t\right)$ which is $\overline{f}\left(s\right)$. In the case of damping oscillation, $\overline{f}\left(s\right)$ takes the form ${\displaystyle \overline{f}\left(s\right)=\frac{\omega }{{(s+\omega )}^2+{\omega }^2}}$. In order to plot the force with respect to time, we seek the inversion of the expression in Equation (37) numerically.

Case 2: Sine oscillation

In this case, we assume that the porous sphere oscillates with the velocity $U\left(t\right)=U_{0}f\left(t\right)$, but with $f\left(t\right)=sin\omega t$, where $\omega$ is the frequency of oscillations. As mentioned in case 1, the drag force is obtained in Laplace domain and it depends on the coefficients coefficients $A_1$ and $A_4$ and the Laplace transform of $f\left(t\right)$ which is $\overline{f}\left(s\right)$. In the sine oscillation case, $\overline{f}\left(s\right)$ takes the form ${\displaystyle \overline{f}\left(s\right)=\frac{\omega }{s^2+{\omega }^2}}$.

Case 3: Impulsive motion

In this case, the porous sphere starts to translate suddenly with the velocity $U\left(t\right)=U_{0}H\left(t\right)$, where $H\left(t\right)$ is Heaviside step function or unit step function. The unit step function is defined as

$$H\left(t\right)=\left\{\begin{array}{cc}0& t<0\\ 1& t\ge 0\end{array}\right.$$

Note that the coefficients $A_1$ and $A_4$ appearing in the force expression obtained in Equation (37) are in Laplace domain and they contain Laplace transform of unit step function $H\left(t\right)$.

5. Numerical Results and Discussions

Here, we present some graphical representations of the streamlines for damping sine oscillation and unit step function.

Figure 1 presents streamline patterns for damping oscillations for various values of fractional parameter $\left(\alpha \right)$ at fixed $T_u=0.5,\beta =0.6,Da=0.3,\omega =1$ and $t=0.5$. It can be seen that at smaller values of $\alpha$, closed streamlines are formed inside and outside of the porous sphere (see Figure 1 at $\alpha =0.01$ and $\alpha =0.1$). At $\alpha =0.1$, the streamlines in the vicinity of the porous sphere are clustered and velocity is highest around the porous sphere. This is observed because the fluid in the vicinity of the porous sphere gets accelerated due to motion of the porous sphere. The velocity decays as we move away from the porous sphere in the y-direction, because far from the porous sphere the magnitude of disturbance to the ambient fluid reduces at smaller value of $\alpha$. Values of $\alpha <0.5$ do not affect the streamlines much. At $\alpha =0.5$, the streamlines near the porous sphere get deformed in the x-direction. In addition, the velocity is maximum in the vicinity of the porous sphere and also far from porous sphere where the streamlines are clustered. At the values $\alpha \ge 0.5$, the effect of unsteadiness is more profound as it contributes to the momentum leading to an increase in velocity far from the porous sphere. It is interesting to note that at $\alpha =1$, closed streamlines are formed in the vicinity of the porous sphere along the x-axis. In addition, the streamlines inside the porous sphere close to the center form parabolic shape which is symmetric about the x-direction. In fact, as $\alpha$ grows the disturbance far from the center of the sphere grows.

Figure 2 presents streamline patterns for unit step function for various values of fractional parameter $\left(\alpha \right)$ at fixed $T_u=0.5,\beta =0.6,Da=0.3,t=0.5$ and $\omega =1$. It may be noticed that at the smaller value of $\alpha =0.01$, the streamlines are undeformed and clustered in the vicinity of the porous sphere leading to maximum velocity on the surface which decays far from the sphere surface. As seen in Figure 2, the smaller values of $\alpha <0.5$ do not affect the streamlines much. At $\alpha \ge 0.5$, streamlines start getting to deform in the vicinity of the porous sphere and also inside the porous sphere. It may be noted that in case of unit step function, the magnitude of disturbance far from the center of the sphere grows as $\alpha$ grows, but at $\alpha =1$, more fluid seeps through the porous sphere forming streamlines inside it as compared to the damping sine oscillation at $\alpha =1$ (see Figure 2 at $\alpha =1$).

Figure 3, Figure 4, Figure 5 and Figure 6 present the non-dimensional drag force, ${\tilde{F}}_z^{\alpha }$, for the three cases mentioned above for different values of the fractional order parameter $\left(\alpha \right)$, Darcy number $\left(Da\right)$, unsteadiness parameter $\left(T_u\right)$ and the stress jump coefficient $\left(\beta \right)$.

Figure 3a–c demonstrate variations of the dimensionless drag force exerted on the spherical surface for different values of the fractional derivative parameter $\alpha$ by fixing the other parameters such as Darcy number $\left(Da\right)$, unsteadiness parameter $\left(T_u\right)$ and the stress jump coefficient $\left(\beta \right)$. In Figure 3a–c, $\beta ,T_u$ and $Da$ are kept at 0.6, 0.5 and 0.01, respectively. It can be observed that the fractional order parameter $\left(\alpha \right)$ makes an effective contribution to the drag. An increase in this parameter increases the values of the drag force. The classical case of non-fractional derivative is achieved when the parameter $\alpha$ is assigned the value of one. It may be concluded that use of fractional order derivative in the Navier–Stokes equation leads to the reduction of drag force on spherical objects. Physically, the smaller the values of $\alpha$, the more rapidly the shear decays which results in reduction of drag force. This may be due to the fractional order derivatives which examines the complete description of the memory effectively. It may also be observed that in the case of damping sine oscillations, drag force tends to zero as $\omega t\to \infty$ which is large time behavior. In this case, the porous sphere does not experience any force arising due to damping effect. It is interesting to see that the particle experiences a reverse force at relatively smaller values of $\omega t$ (Figure 3a). In case of sine oscillation, the force profile exhibits oscillatory behavior and it never approaches zero (Figure 3b). Unlike the damping sine and sine oscillations, the porous sphere does not experience negative force in case of unit step function.

Figure 4a–c exhibit variations of drag force versus the time parameter for different values of stress jump coefficient $\left(\beta \right)$ at fixed $\alpha =0.5$, $T_u=0.5$ and $Da=0.01$. It can be noticed that an increase in this coefficient results in a decrease in the drag force. This may be due to the fact that increase in the stress jump coefficient $\left(\beta \right)$ forces more fluid to seep through the porous sphere before it slips from the surface which leads to a reduction in the drag force. Figure 5a–c depict the variation of drag force versus the time parameter for different values of unsteadiness parameter $\left(T_u\right)$ at fixed $\alpha =0.5$, $\beta =0.6$ and $Da=0.01$. It may be noted that increase in unsteadiness parameter contributes to the momentum of the fluid flow which further helps the shear stress acting on the porous sphere leading to an increase in the drag force experienced by the porous sphere.

Figure 6a–c present variation of drag force versus the time parameter for different values of Darcy number $\left(Da\right)$ at fixed $\alpha =0.5$, $\beta =0.6$ and $T_u=0.5$. It can readily be seen that an increase in Darcy number $\left(Da\right)$ reduces the drag force. This occurs due to the fact that increase in Darcy number increases the permeability of the porous medium which helps ease the fluid flow inside the porous particle and this leads to a decrease in the drag force. It is to be noted that the force reflects the pattern of the basic flow, for example if the basic flow is damping sine, the force exhibits damping behavior and eventually becomes zero at larger time. If the basic flow is oscillatory, the force exhibits oscillatory behavior. In case of the basic flow taken as unit step function, the force obtained is large in magnitude at $\omega t=0$ and the magnitude decreases as $\omega t$ grows. It can be concluded that the use of fractional order derivatives exhibits the history (‘memory’) of the basic flow in the obtained physical quantities of interest.

6. Conclusions

This paper presents a semi analytical solution for the problem of unsteady translational motion of a porous sphere in an incompressible viscous fluid. The fluid flow exterior and interior of the porous sphere is governed by the modified Navier–Stokes equation with fractional time derivative and the modified Brinkman equation with fractional time derivative. The Laplace transform technique is used and the inversion is obtained with the aid of numerical inversion formula of Laplace transform. Streamline patterns for damping sine and unit step function are presented. It was observed that increase in the fractional parameter $\alpha$ increases the disturbance far from the center of the porous sphere in both the cases. However, it was noticed that at smaller values of $\alpha$, the streamlines do not deform outside and inside the porous sphere. It may also be noted that in case of unit step function, at $\alpha =1$, more fluid seeps through the porous sphere resulting in formation of clustered streamlines leading to increase in velocity inside the porous sphere compared to the case of damping sine oscillation. The drag force formula is obtained analytically in Laplace domain and is inverted numerically. Variation of the drag force with various parameters such as fractional parameter $\left(\alpha \right)$, stress jump coefficient $\left(\beta \right)$, unsteadiness parameter $\left(T_u\right)$ and Darcy number $\left(Da\right)$ is discussed graphically. It is concluded from the numerical results that the fractional order parameter, $\alpha$, contributes considerably on the drag force. It is observed that the values of the drag force increase with increase of this parameter. In addition, it is observed that drag force decreases with stress jump coefficient and Darcy number whereas it increases with unsteadiness parameter.