Transient pressure driven flow in an annulus partially filled with porous material: Azimuthal pressure gradient

Consider the motion of transient, fully developed and incompressible circumferential flow due to Azimuthal pressure gradient between two co-axial-horizontal cylinders containing fluid and porous layer separated by a permeable thin interface with fluid occupying the interval $r_{i} \leq r^{\prime} \leq d^{\prime}$ while the interval $d^{\prime} \leq r^{\prime} \leq r_{0}$ is occupied by a fluid-saturated porous material of uniform permeability.

 and  are the radii of the inner and the outer cylinder respectively as shown in Fig. 1. The circumferential velocity component of the fluid is proved to be a function of the radial coordinate  and time  only. Following the work of Tsangaris and Vlachakis [11], the fundamental model governing the present physical situation in dimensional form can be written as follows.

Figure 1. Flow configuration and coordinate system

$\frac{\partial u_{p}^{\prime}}{\partial t^{\prime}}=v_{e f f}\left[\frac{\partial^{2} u_{p}^{\prime}}{\partial r^{\prime 2}}+\frac{1}{r^{\prime}} \frac{\partial u_{p}^{\prime}}{\partial r^{\prime}}-\frac{u_{p}^{\prime}}{r^{\prime 2}}\right]-\frac{v}{k^{\prime}} u_{p}^{\prime}-\frac{1}{\rho} \frac{\partial p}{\partial \theta} \frac{1}{r^{\prime}}$ for $r_{i} \leq r^{\prime} \leq d^{\prime}$   (1)

$\frac{\partial u_{f}^{\prime}}{\partial t^{\prime}}=v\left[\frac{\partial^{2} u_{f}^{\prime}}{\partial r^{\prime 2}}+\frac{1}{r^{\prime}} \frac{\partial u_{f}^{\prime}}{\partial r^{\prime}}-\frac{u_{f}^{\prime}}{r^{\prime 2}}\right]-\frac{1}{\rho} \frac{\partial p}{\partial \theta} \frac{1}{r^{\prime}} \quad$ for $\quad d^{\prime} \leq r^{\prime} \leq r_{0}$   (2)

The initial and the no-slip condition at the surfaces for the problem under consideration are

$t^{\prime} \leq 0 : u_{p}^{\prime}=u_{f}^{\prime}=0$ for $r_{i} \leq r^{\prime} \leq r_{0}$

$t^{\prime}>0 :\left[\begin{array}{l}{u_{p}^{\prime}=0 \text { at } r^{\prime}=r_{i}} \\ {u_{f}^{\prime}=0 \text { at } r^{\prime}=r_{0}}\end{array}\right]$   (3)

With the dimensional matching condition at the interface given as

$t^{\prime}>0 :\left[\begin{array}{l}{u_{p}^{\prime}=u_{f}^{\prime}=u_{i}^{\prime}} \\ {\text {veff}\left[\frac{\partial u_{p}^{\prime}}{\partial r^{\prime}}-\frac{u_{p}^{\prime}}{r^{\prime}}\right]-v\left[\frac{\partial u_{f}^{\prime}}{\partial r^{\prime}}-\frac{u_{f}^{\prime}}{r^{\prime}}\right]=\frac{\beta v}{\sqrt{k^{\prime}}} u_{p}^{\prime}}\end{array}\right]$ at $r^{\prime}=d^{\prime}$

Introducing the following dimensionless quantities in Eqs. (1) – (3)

$R=\frac{r^{\prime}}{r_{i}} ; t=\frac{v t^{\prime}}{r_{i}^{2}} ; \gamma=\frac{v_{e f f}}{v} ; D a=\frac{r_{i}^{2}}{k^{\prime}} ; U_{p}=U_{f}=\frac{\left(u_{p}^{\prime} ; u_{f}^{\prime}\right)}{U_{0}}$

$U_{0}=-r_{i} \frac{\partial p}{\partial \theta} \frac{1}{\rho v} ; d=\frac{d^{\prime}}{r_{i}} ; \lambda=\frac{r_{0}}{r_{i}} ; U_{i}=\frac{u_{i}^{\prime}}{U_{0}}$   (4)

Equations (1) and (2) can be written in dimensionless form as

$\frac{\partial U_{p}}{\partial t}=\gamma\left[\frac{\partial^{2} U_{p}}{\partial R^{2}}+\frac{1}{R} \frac{\partial U_{p}}{\partial R}-\frac{U_{p}}{R^{2}}\right]-\frac{U_{p}}{D a}+\frac{1}{R}$ for $1 \leq R \leq d$   (5)

$\frac{\partial U_{f}}{\partial t}=\left[\frac{\partial^{2} U_{f}}{\partial R^{2}}+\frac{1}{R} \frac{\partial U_{f}}{\partial R}-\frac{U_{f}}{R^{2}}\right]+\frac{1}{R} \quad$ for $\quad d \leq R \leq \lambda$   (6)

Subject to the following non-dimensional initial and boundary condition

$t \leq 0 : U_{p}=U_{f}=0$ for $1 \leq R \leq \lambda$

$t>0 :\left[\begin{array}{l}{U_{p}=0 \text { at } R=1} \\ {U_{f}=0 \text { at } R=\lambda}\end{array}\right]$   (7)

With the non-dimensional matching condition at the interface given as

$t>0 :\left[\begin{array}{l}{U_{p}=U_{f}=U_{i}} \\ {\gamma\left[\frac{\partial U_{p}}{\partial R}-\frac{U_{p}}{R}\right]-\left[\frac{\partial U_{f}}{\partial R}-\frac{U_{f}}{R}\right]=\frac{\beta}{\sqrt{D a}} U_{p}}\end{array}\right]$ at $R=d$ 

The expression of equation $(5)$ to $(7)$ in Laplace domain using the Laplace transform technique $\overline{U}_{p}(\mathrm{R}, \eta)=\int_{0}^{\infty} U_{p}(R, t) e^{-\eta t} d t$ and $\overline{U}_{f}(R, \eta)=\int_{0}^{\infty} U_{f}(R, t) e^{-\eta t} d t$ (Where $\eta>0$ is the Laplace parameter) subject to initial condition yield:

$\frac{d^{2} \overline{U}_{p}}{d R^{2}}+\frac{1}{R} \frac{d \overline{U}_{p}}{d R}-\frac{\overline{U}_{p}}{R^{2}}-\frac{1}{\gamma}\left[\frac{1}{D a}+\gamma\right] \overline{U}_{p}=-\frac{1}{\gamma \eta R}$   (8)

$\frac{d^{2} \overline{U}_{f}}{d R^{2}}+\frac{1}{R} \frac{d \overline{U}_{f}}{d R}-\frac{\overline{U}_{f}}{R^{2}}-\eta \overline{U}_{f}=-\frac{1}{\eta R}$   (9)

Subject to the following non-dimensional boundary condition

$t>0 :\left[\begin{array}{l}{\overline{U}_{p}=0 \text { at } R=1} \\ {\overline{U}_{f}=0 \text { at } R=\lambda}\end{array}\right]$   (10)

With the non-dimensional matching condition at the interface given as

$t>0 :\left[\begin{array}{c}{\overline{U}_{p}=\overline{U}_{f}=U_{i}} \\ {\gamma\left[\frac{d \overline{U}_{p}}{d R}-\frac{\overline{u}_{p}}{R}\right]-\left[\frac{d \overline{U}_{f}}{d R}-\frac{\overline{U}_{f}}{R}\right]=\frac{\beta}{\sqrt{D a}} \overline{U}_{p}}\end{array}\right]$ at $R=d$ 

Equations (8) and (9) can be reduced by using the following transformations respectively

$\overline{U}_{p}=\overline{U}_{p h}+\frac{1}{\eta\left(\frac{1}{D a}+\eta\right)_{R}} \quad$ and $\overline{U}_{f}=\overline{U}_{f h}+\frac{1}{\eta^{2} R}$   (11)

applying equation (11) on equations (8) and (9), the Laplace domain solution in terms of the modified Bessel function obtained are given by;

$\overline{U}_{p}=C_{1} I_{1}(\delta R)+C_{2} K_{1}(\delta R)+\frac{1}{\eta\left(\frac{1}{D a}+\eta\right)_{R}}$   (12)

$\overline{U}_{f}=C_{3} I_{1}(\sqrt{\eta} R)+C_{4} K_{1}(\sqrt{\eta} R)+\frac{1}{\eta^{2} R}$   (13)

where $\delta=\sqrt{\frac{1}{\gamma}\left(\frac{1}{D a}+\eta\right)}$ and $I_{1}, K_{1}$ are the modified Bessel function of first and second kind respectively of order $1 .$

Applying the boundary condition (10) on equation (12) and (13) the constants C₁, C₂, C₃, C₄ and the expression for the interfacial velocity U_i are obtained as

$C_{1}=x_{2}-U_{i} x_{3}, C_{2}=U_{i} x_{4}-x_{5}, C_{3}=x_{7}-U_{i} x_{8}, C_{4}=U_{i} x_{10}-x_{9}$

$U_{i}=\frac{x_{13}+x_{12} x_{5}-x_{11} x_{2}+x_{14} x_{7}-x_{15} x_{9}-x_{16}}{x_{14} x_{8}+x_{15} x_{10}-x_{11} x_{3}-x_{12} x_{4}-\frac{\beta}{\sqrt{D a}}}$   (14)

The constant $x_{1}, x_{2}, x_{3}, x_{4}, \ldots x_{16}$ appearing in the above equations are defined in the Appendix

2.1 Skin friction

The skin friction at the outer surface of the inner cylinder (R=1) and at the inner surface of the outer cylinder R=λ in Laplace domain is given respectively by;

$\overline{\tau}_{1}=R \frac{d}{d R}\left.\left(\frac{\overline{U}_{p}}{R}\right)\right|_{R=1}=\delta\left[C_{1} I_{2}(\delta)-C_{2} K_{2}(\delta)-\frac{2}{\delta \eta\left(\frac{1}{D a}+\eta\right)}\right]$   (15)

$\overline{\tau}_{\lambda}=-R \frac{d}{d R}\left(\frac{\overline{U}_{f}}{R}\right)| |_{R=\lambda}=\sqrt{\eta}\left[C_{4} K_{2}(\sqrt{\eta} \lambda)-C_{3} I_{2}(\sqrt{\eta} \lambda)+\right.\frac{2}{\eta \lambda^{2}}$    (16)

Equations (12) - (16) are then inverted to determine their solutions in the time domain. However, literature survey conducted indicated that there are no simple ways to invert Laplace domain solutions in modified Bessel function to the time domain. Hence, we adopt a numerical procedure used in the work of Jha and Yusuf [31], Jha and Odengle [16], Jha and Apere [32] as well as that of Khadrawi and Al-Nimr [33] which were founded on the Riemann-sum approximation. According to this method, any function in the Laplace domain can be inverted to the time domain as follow:

$Z(R, t)=\frac{e^{\varepsilon t}}{t}\left[\frac{1}{2} \overline{Z}(R, \varepsilon)+\operatorname{Re} \sum_{k=1}^{M} \overline{Z}(R, \varepsilon+\right.\frac{i k \pi}{t} )(-1)^{k} ], 1 \leq R \leq \lambda$   (17)

where Re represents the real part of the term with summation and $i=\sqrt{-1} . \mathrm{M}$ is the number of terms used in the Riemann-sum approximation and

 is the real part of the Bromwich contour that is used in inverting Laplace transforms. The Riemann-sum approximation for the Laplace inversion involves a single summation for the numerical process its accuracy generally depends on the value of ε and M. Following the pioneering work of Tzou [34], the value of   that best satisfied the result is 4.7.

2.2 Validation of the method

To validate the correctness of the Riemann-sum approximation scheme adopted in inverting Equation (12)–(16), we find analytical solution of the steady state circumferential velocity, which should coincide with the transient solution at large time. This can be obtained by setting $\frac{\partial( )}{\partial t}=0$ in equations $(2.5)$ and $(2.6)$ to obtain the following dimensionless ordinary differential equations

$\frac{d^{2} U_{p}}{d R^{2}}+\frac{1}{R} \frac{d U_{p}}{d R}-\frac{U_{p}}{R^{2}}-\frac{U_{p}}{\gamma D a}=-\frac{1}{\gamma R}$   (18)

$\frac{d^{2} U_{f}}{d R^{2}}+\frac{1}{R} \frac{d U_{f}}{\mathrm{d} R}-\frac{U_{f}}{R^{2}}=-\frac{1}{R}$   (19)

Subject to the boundary conditions

$U_{p}=0$ at $R=1$   

$U_{f}=0$ at $R=\lambda$   (20)

With the non-dimensional matching condition at the interface given as

$\left[\begin{array}{l}{U_{p}=U_{f}=U_{i}} \\ {\gamma\left[\frac{d \nu_{p}}{d R}-\frac{U_{p}}{R}\right]-\left[\frac{d U_{f}}{d R}-\frac{U_{f}}{R}\right]=\frac{\beta}{\sqrt{D a}} U_{p}}\end{array}\right]$ at $R=d$   (21)

The solution of equation $(18)$ and $(19)$ are given by applying appropriate $\quad$ transformation $U_{p}=U_{p h}+\frac{D a}{R} \quad$ and $\quad R=e^{\xi}$ respectively, to yield

$U_{p}=C_{5} I_{1}(\zeta R)+C_{6} K_{1}(\zeta R)+\frac{D a}{R}$   (22)

$U_{p}=C_{7} R+\frac{c_{8}}{R}-\frac{R l n(R)}{2}$   (23)

where $\zeta=\frac{1}{\sqrt{\gamma D a}}$

 Applying the boundary condition $(20)$ and the condition at the interface $(21)$ on equations $(22)$ and $(23),$ the unknown constants $C_{5}, C_{6}, C_{7}, C_{8}$ as well as the interfacial velocity $U_{i}$ are given as

$C_{5}=y_{4}-U_{i} y_{5}, C_{6}=y_{6}+U_{i} y_{7}, C_{7}=y_{11}+U_{i} y_{12}, C_{8}=y_{8}+U_{i} y_{9}$

$U_{i}=\frac{\frac{2 y_{8}}{d^{2}}+y_{1} y_{4}-y_{10}-y_{2} y_{6}}{y_{1} y_{5}+y_{2} y_{7}-\frac{2 y_{9}}{d^{2}}}$   (24)

The constant $y_{1}, y_{2}, y_{3}, y_{4}, \ldots y_{12}$ appearing in the above equations are defined in the Appendix

In an attempt to further establish the accuracy of the numerical procedure used in this research work, we used the implicit finite difference method to solve the dimensionless Equations (5) and (6) with the initial and the no slip boundary conditions (7). Using this numerical scheme, an excellent agreement was found in comparison with the values obtained from Riemann-sum approximation approach at large time as well as the exact solution of the steady state velocity. It is good to note that, the numerical values obtained using the implicit finite difference method at transient state also agrees considerably with the ones obtained using the Riemann-sum approximation approach for a small value of time. (See Table 1 & 2).

Table 1. Numerical value for the steady-state velocity obtained using the Riemann-sum approximation, implicit finite difference and exact solution (γ=1, β=0, λ=2, d=1.5)

R	Da	Riemann-sum approximation	Implicit finite difference	Exact solution
1.2	0.01	0.0078	0.0078	0.0078
	0.10	0.0370	0.0370	0.0370
	1.00	0.0566	0.0566	0.0566
1.4	0.01	0.0116	0.0116	0.0116
	0.10	0.0514	0.0514	0.0514
	1.00	0.0769	0.0768	0.0769
1.6	0.01	0.0265	0.0266	0.0265
	0.10	0.0544	0.0543	0.0543
	1.00	0.0715	0.0715	0.0715
1.8	0.01	0.0235	0.0235	0.0235
	0.10	0.0366	0.0366	0.0366
	1.00	0.0447	0.0446	0.0446

The semi-analytical solution of the time dependent flow in an annulus partially filled with porous material due to sudden application of constant circumferential (Azimuthal) pressure gradient  is considered. The annular gap is filled with a clear fluid and a fluid-saturated porous material of uniform permeability separated by a thin porous layer, with the clear fluid occuping the region near the inner surface of the outer cylinder. Laplace transform technique is used to solve the non-dimensional transport equations. The solution in Laplace domain is then inverted to time-domain using Riemann-sum approximation approach. In order to establish the accurary of the Riemann-sum approximation approach, a number of tables have been presented which shows that at large time, Riemann-sum approximation approach largely agrees with the steady state solution obtained exactly. To further demostrate reliability of this approach, a numerical scheme – the implicit finite defference is also used to solve the governing equations. The numerical values obtained from this two approaches (See table 1-4) indicates that they both match considerably at transient and steady state. The influence of a number of selected controlling parameters entering into the model is discused with the aid of line graphs. The main conclusions of the research work are

1. Generally, the unified velocity  is an increasing function of time (t) till it attains steady state velocity.

2. A thin clear fluid region boosts τ₁but has no influence on τ_λ. While thin porous region enhances skin friction at both surfaces.

3. There is a clear indication that porous region is more permeable with greater fluid motion arising in the porous region for a large value of Darcy number, this transitively enhanced the skin friction at the surfaces, The reverse trend occurs for a small value of Darcy number.

4. Positive value of β is observed to enhance interfacial velocity U_i, while large value of d (as d tends to the radii ratio) retards U_i. Furthermore, a linear relationship (See Fig. 8) is seen at steady state with negative gradient.