On Modeling Laminar Flow through Variable Permeability Transition Layer

1. Introduction

Driven by a host of well-documented applications in natural and industrial settings, including movement of water, oil and gas in ground layers, lubrication theory, blood flow in biological tissues, and various other industrial applications (cf. Bejan and Nield [1] and the references therein) fluid flow through and over porous layers continues to receive unparalleled attention in the porous media literature. Depending on the porous microstructure of the medium, the presence or absence of macroscopic, solid boundaries on which the no-slip condition is imposed, and the importance of inertial effects in the medium, the flow of an incompressible fluid is governed by the continuity equation

$$\nabla ·\overrightarrow{v}=0$$

(1)

and a general momentum equation which, in the absence of body force, is written in the following form when permeability is constant [2,3]

$$\rho {\overrightarrow{v}}_t+[\beta \left(\delta -1\right)+1](\overrightarrow{v}·\nabla )\overrightarrow{v}]=-\nabla p+\mu [\beta \left(\vartheta -1\right)+1]{\nabla }^2\overrightarrow{v}-\beta \left\{{\displaystyle \frac{\mu }{k}}\overrightarrow{v}+{\displaystyle \frac{\rho \sigma C_d}{\sqrt{k}}}\overrightarrow{v}\left|\overrightarrow{v}\right|\right\}$$

(2)

where $\overrightarrow{v}$ is the fluid velocity vector field, $k$ is the permeability tensor, $p$ is the pressure, $\rho$ is the fluid density, $\mu$ is the base fluid viscosity, $\vartheta ={\displaystyle \frac{{\mu }_e}{\mu }}$ is the relative viscosity, ${\mu }_e$ is the effective viscosity (viscosity of the fluid in the porous medium), $\delta$ is the convective acceleration parameter, $\beta$ and $\sigma$ are binary parameters, and $C_d$ is the Forchheimer drag coefficient.

The parameters $\beta ,\delta ,\sigma$ take on the values $0$ or $1$. By assigning the values $0$ or $1$ to parameters $\beta ,\delta ,\sigma$, the main models of fluid flow through porous media and in free-space are recovered. For instance, if $\beta =0$, the flow is in free-space and equation (2) reduces to the Navier-Stokes equations. If $\beta =1$, the flow is in the porous medium, and the following flow models (Table 1) are recovered, [3].

Flow over porous layers involves flow through a free-space channel over a porous layer, or flow through a composite of porous layers of different porosities and permeabilities (cf. [4,5,6,7,8] and the references therein). Its direct influence on heat and mass transfer processes in various physical settings mandates accurate modelling of the flow processes and the associated interfacial conditions between flow domains, [9,10,11,12,13,14]. Based on equation (1), uni-directional flow in channels and porous layers is governed by the following equation, wherein $u$ is the tangential velocity component:

$$-p_x+\mu [\beta \left(\vartheta -1\right)+1]u_{yy}-\beta \left\{{\displaystyle \frac{\mu }{k}}u+{\displaystyle \frac{\rho \sigma C_d}{\sqrt{k}}}{u}^2\right\}$$

(3)

Equation (3) provides a single-equation for simulation of flow through and over porous layers by simply switching parameters.

Experimental and theoretical studies in the field of flow over porous layers span more than six decades, and include the development of an expansive body of knowledge that initially focused on the use of Darcy’s equation as a model of flow through the porous layer and resulted in the development of Beavers and Joseph slip hypothesis, [15], and the following condition at the interface, $(y=0)$, between a free-space layer and the porous layer

$$u_y\left(x,0^{+}\right)={\displaystyle \frac{\alpha }{\sqrt{k}}}[u\left(x,0^{+}\right)-u_D]$$

(4)

where $u_D$ is the Darcy velocity, and $\alpha$ is an empirical slip parameter with a range of values 0.1 to 4, (cf. [1,15,16,17].

The Beavers and Joseph condition continues to be addressed by various authors (cf. [1,18,19,20,21] and the references therein). Concerns about the low order of Darcy’s equation and the slip parameter in Beavers and Joseph’s experiments resulted in the development of interfacial conditions centred on velocity and shear stress continuities, and the use of Brinkman’s equation, [22], as a viable alternative to Darcy’s equation. Furthermore, if Darcy’s equation is replaced by another flow model, such as the Forchheimer equation, then the slip parameter $\alpha$ needs to be adjusted, as discussed by Roach and Hamdan, [23]. Neale and Nader, [19], provided the following conditions at the assumingly sharp interface between a free-space channel$0<y<h)$ and a Brinkman porous layer (-h<y<0), and identified the slip parameter $\alpha$ with $\sqrt{\vartheta }$:

$$u\left(x,0^{+}\right)=u\left(x,0^{-}\right)$$

(5)

$$u_y\left(x,0^{+}\right)={\vartheta u}_y\left(x,0^{-}\right)$$

(6)

Brinkman’s equation, [22], continues to receive attention by various authors, whether to discuss its validity or its limitations (cf. [1,18,24,25,26] and the references therein).

It is generally agreed that, in the flow over a porous layer, if the porous layer is of infinite depth, then the flow in the porous layer is governed by Darcy’s equation, while if the porous layer is of finite depth, then Brinkman’s equation should be used (cf. Rudraiah [7] and the references therein).

Ochoa and Whitaker [20,21] developed a stress jump condition at the interface, which takes into account momentum transfer at the interface, and is of the form

$${{\mu }_{e}u}_y\left(x,0^{+}\right)-{\mu u}_y\left(x,0^{-}\right)=-{\displaystyle \frac{\gamma }{\sqrt{k}}}u$$

(7)

where $\gamma$ is an adjustable jump parameter. The absence of definite values $\gamma$ hinders the efficient use of the Ochoa-Whitaker jump stress condition , hinders its accurate use, [28].

Many of the interfacial conditions and models in use assume the existence of a sharp interface, and do not account for permeability discontinuity at the interface between free-space and the porous layer, or between two porous layers of different permeabilities. Extensive analyses have been carried out to point the need for interfacial modeling that takes into account gradual changes in permeability from its finite value in the porous medium to its approach to infinity in the absence of a porous matrix, or a gradual change in permeability from one finite value to another different finite value in the case of flow through a porous layer underlain by another porous layer, [29].

Recent advancements in the field, however, concluded that permeability continuity is better handled through the introduction of a thin transition region between a homogeneous porous layer and free-space, or between two homogeneous porous layers. The concept of the transition layer has been implemented and validated in the analysis provided by number of authors [9,13,14,18,28,30,31,32,33,34]. The transition layer is characterized by its variable permeability that smoothly changes between its value in the homogeneous porous sediment and infinity in the free-space channel, or between two finite permeability values in two adjacent homogeneous porous layers. The crux of the problem is then narrowed down to modeling flow through a variable permeability porous layer the flow through which is governed by Brinkman’s equation. This has given an impetus to Brinkman’s equation and its usefulness.

It is worth noting that the use of the transition layer is not restricted to flow over a porous layer. Rather, three transition layer scenarios are possible in this type of tri-layer analysis. The first is a transition layer between a Navier-Stokes channel and a constant permeability porous layer of infinite depth. In this case, the flow in the constant permeability layer is governed by Darcy’s law.

The second is a transition layer between the Navier-Stokes channel and a homogeneous porous layer of finite depth. In this case, the flow in the constant permeability finite layer is governed by Brinkman’s equation. The third involves a transition layer inserted between two Darcy layers, or between two homogeneous Brinkman layers, or a combination of both. In all of these cases, the mathematical problem reduces to that of a Poiseuelle-type flow configuration, or a uni-directional flow problem in a tri-layer flow domain.

Many investigations and models exist for flow through variable permeability porous layers, (cf. [35,36,37,38] and the references therein). Roach and Hamdan, [39], provided a review and a classification of the available variable permeability models.

A most suitable model was introduced in the work or Nield and Kuznetsov [30], where the reciprocal of the permeability varies across the transition layer according to $K=a{k}_{0}H/y$, where $H$ is the overall tri-layer thickness, $a$ is related to the thickness of the transition layer, $K$ is the transition layer permeability, and $k_0$ is the constant permeability of the underlying porous layer. This choice of permeability distribution reduces Brinkman’s equation in the transition layer to an Airy’s inhomogeneous ordinary differential equation of the form:

$$u_{yy}-{\displaystyle \frac{1}{\vartheta Da}}\left({\displaystyle \frac{y-\xi }{\zeta }}\right)u+{\displaystyle \frac{1}{\vartheta }}=0;\xi <y<\eta ;\zeta =\eta -\xi$$

(8)

where $\zeta$ is the transition layer thickness, $u$ is the flow velocity in the transition layer, and $Da={\displaystyle \frac{k_0}{H^2}}$ is the Darcy number that is defined in terms of a constant reference permeability, $k_0$. Equation (8) is the well-known Airy’s inhomogeneous equation.

In the process of their solution, Nield and Kuznetsov [30] introduced a new integral function, $N_i\left(y\right)$, to express the particular solution to the resulting Airy’s inhomogeneous equation. The function $N_i\left(y\right)$ is referred to in the literature as the Nield-Kuznetsov function, [40], and is defined in terms of Airy’s functions of the first- and second-kind, $A_i\left(y\right)$ and $B_i\left(y\right)$, respectively, as

$$N_i\left(y\right)=A_i\left(y\right){\int }_0^y{B}_i\left(t\right)dt-B_i\left(y\right){\int }_0^y{A}_i\left(t\right)dt$$

(9)

General solution to (8) is thus expressed as:

$$u_2\left(y\right)=c_2{A}_i\left(\lambda (y-\xi )\right)+d_2{B}_i\left(\lambda (y-\xi )\right)+{\displaystyle \frac{\pi }{\vartheta {\lambda }^2}}{N}_i\left(\lambda \right(y-\xi \left)\right)$$

(10)

where $\lambda ={\displaystyle \frac{1}{\sqrt[3]{\vartheta Da\left(\eta -\xi \right)}}}$The Nield-Kuznetsov $N_i\left(y\right)$ function and its role in the transition layer was successfully implemented by Tao et al. [28], who presented elegant and detailed analysis of the effects of the transition layer parameters involved in the modelling process, where they employed two specific permeability functions to express permeability variations. Tao et al. [28] provided details of the amplifications in the transition layer permeability and emphasized the need for an exact description of the permeability in the transition layer and a precise permeability forecasting function to model the gradual change of permeability along the transition layer region to better match the velocity distribution with experimental results.

Tao et al. [28] employed ascending series representations in the evaluations of ${A_i\left(y\right),B_i\left(y\right)\mathrm{a}\mathrm{n}\mathrm{d}N}_i\left(y\right)$. Nield and Kuznetsov [30] evaluated then using asymptotic series expansions and reported that the smallest value of $Da$ they could compute their results for was $Da=0.0002$, and stated that as $Da$ tends to zero, $\lambda ={\displaystyle \frac{1}{\sqrt[3]{\vartheta Da\left(\eta -\xi \right)}}}$ becomes large and results in singularity. While there is no claim that we know the reason as to why this happened, it is quite possible that the above-mentioned use of asymptotic series (valid for large values of $y$) as opposed to ascending series (valid for all values of $y$) played a role. We also point out that in the expression for $\lambda$, the transition layer thickness, $\eta -\xi$, is a positive, finite quantity less than unity, and that $\vartheta ={\displaystyle \frac{{\mu }_e}{\mu }}$ is also finite. If we write the coefficient of $u$ in equation (8) as ${\displaystyle \frac{y-\xi }{\vartheta Da(\eta -\xi )}}$, then as $y\to \eta$, and $Da\to 0$, the expression ${\displaystyle \frac{y-\xi }{\vartheta Da(\eta -\xi )}}$ becomes excessively large, unless the numerator in the expression approaches zero at the same rate as $Da$ approaches zero. One possibility is to introduce a permeability model that leads to ${(y-\xi )}^n$ in the numerator, where $n$ is a positive integer. By adjusting the value of $n$, it might be possible to obtain solutions for small values of $Da$. In the current work, the smallest value of Darcy number that we could use in computations involving a thin transition layer and the use of ${(y-\xi )}^n$ was $Da=0.0000005$. Furthermore, the suggested choice of permeability function might prove effective in reducing permeability amplifications in the transition layer that Tao et al. [28] discussed.

This observation, together with the success attained by Nield and Kuznetsov, [30], and Tao et al. [28], in advancing the state of knowledge of flow in the transition layer, motivates the introduction of other integral functions to model variations of permeability in the transition layer. In fact, in a recent study where a porous layer was embedded between two fluid layers, Abu Zaytoon and Hamdan, [41], employed a permeability function that reduced Brinkman’s equation to Weber’s inhomogeneous equation, the solution to which they obtained in terms of the parametric Nield-Kuznetsov function.

In the current work, we introduce a generalized permeability function, to describe the variable permeability in the transition layer, that offers modeling flexibility and can be reduced to the permeability function used by Nield and Kuznetsov [30]. Brinkman’s equation governing flow in the transition layer can then be reduced to a generalized Airy’s equation (in the sense introduced by Swanson and Headley [42]). An efficient evaluation procedure of the special functions that arise in the solution of the resulting governing equation is presented, and solutions are obtained for a choice of flow parameters. The solution provided by Nield and Kuznetsov [30] is then compared with a special case of the current formulation, which shows agreement between solutions obtained. It is believed that the current approach will provide extended flexibility in modelling permeability variations in the transition zone.

2. Problem Formulation

Consider the flow of an incompressible fluid in the tri-layer model of Figure 1, below, in which layer 1 represents a free-space channel the flow through which is governed by Navier-Stokes equations. Layer 3 is a homogeneous porous layer of constant permeability, the flow through which is governed by Brinkman’s equation. Layer 2 is the transition layer the flow through which is governed by Brinkman’s equation with variable permeability. The flow in the configuration is driven by a common pressure gradient.

Permeability distribution is given as follows, where we take $K_0$ as a constant reference permeability, and assume that flow in the free-space layer 1 has an infinite conductivity to the fluid.

$$\mathrm{Region}1:K_1\to \infty \mathrm{for}\eta H<y^{*}<H$$

(11)

$$\mathrm{Region}2:K_2={\displaystyle \frac{K_0[\left(\eta -\xi \right)H{]}^n}{(\eta H-y^{*}{)}^n}};\xi H<y^{*}<\eta H$$

(12)

$$\mathrm{Region}3:K_3=K_0;0<y^{*}<\xi H$$

(13)

Equations governing the flow in the respective regions are given by:

$$\mu {\displaystyle \frac{d^2{u}_1^{*}}{d{y}^{*2}}}+G=0;\eta H<y^{*}<H$$

(14)

$${\mu }_{e2}{\displaystyle \frac{d^2{u}_2^{*}}{d{y}^{*2}}}-{\displaystyle \frac{\mu }{K_2}}{u}_2+G=0;\xi H<y^{*}<\eta H$$

(15)

$${\mu }_{e3}{\displaystyle \frac{d^2{u}_3^{*}}{d{y}^{*2}}}-{\displaystyle \frac{\mu }{K_3}}{u}_2+G=0;0<y^{*}<\xi H$$

(16)

In the above equations, $G=-{\displaystyle \frac{dp}{dx}}$ is the constant pressure gradient, the quantities $u_i^{*}$ and $K_i$represent the velocity and permeability in region $i$, for $i=\mathrm{1,2},3,$ respectively; $\mu$ is the base fluid viscosity and ${\mu }_{e2}$, ${\mu }_{e3}$ are the effective viscosities, respectively, in regions 2 and 3.

We point out here that in the above formulation, we interchanged layers 1 and 3 in the original formulation by Nield and Kuznetsov, [30]. However, for the sake of comparing results, we solve the flow equations in both configurations.

Following Nield and Kuzentsov [30], we introduce the following dimensionless variables, where $Da={\displaystyle \frac{K_0}{H^2}}$ refers to Darcy number; $y={\displaystyle \frac{y^{*}}{H}}$; $u_i={\displaystyle \frac{\mu }{G{H}^2}}{u}_i^{*}$ for $i=\mathrm{1,2},3$; and $M_i={\displaystyle \frac{{\mu }_{ei}}{\mu }}$ for $i=\mathrm{2,3}$, to transform governing equations (14)-1(6), respectively, into:

$${\displaystyle \frac{d^2{u}_1}{d{y}^2}}+1=0;\eta <y<1$$

(17)

$${\displaystyle \frac{d^2{u}_2}{d{y}^2}}-({\tau }_n{)}^{n+2}{\left(\eta -y\right)}^n{u}_2+{\displaystyle \frac{1}{M_2}}=0;;\xi <y<\eta$$

(18)

$${\displaystyle \frac{d^2{u}_3}{d{y}^2}}-{\displaystyle \frac{1}{M_{3}Da}}{u}_3+{\displaystyle \frac{1}{M_3}}=0;0<y<\xi$$

(19)

where

$${\tau }_n={\displaystyle \frac{1}{\sqrt[n+2]{M_{2}Da(\eta -\xi {)}^n}}}$$

(20)

The dimensionless permeability distributions, ${\displaystyle \frac{K_i}{H^2}},i=\mathrm{2,3}$, can be written in the following dimensionless form in terms of Darcy number, $Da$:

$${\displaystyle \frac{K_2}{H^2}}=Da{\left({\displaystyle \frac{\eta -\xi }{\eta -y}}\right)}^n\mathrm{and}{\displaystyle \frac{K_3}{H^2}}=Da$$

(21)

Equations (17), (18) and (19) must be solved subject to the following boundary and matching conditions

$$u_1\left(1\right)=u_3\left(0\right)=0$$

(22)

$$u_1\left(\eta \right)=u_2\left(\eta \right);u_2\left(\xi \right)=u_3\left(\xi \right)$$

(23)

$${\displaystyle \frac{d{u}_1}{dy}}\left(\eta \right)={\displaystyle \frac{{\mu }_{e2}}{\mu }}{\displaystyle \frac{d{u}_2}{dy}}\left(\eta \right);{\displaystyle \frac{d{u}_2}{dy}}\left(\xi \right)={\displaystyle \frac{{\mu }_{e3}}{{\mu }_{e2}}}{\displaystyle \frac{d{u}_3}{dy}}\left(\xi \right)$$

(24)

We note that if $n=1$, then ${\tau }_n$ corresponds to ${\lambda }_2$ that was used in the work of Nield and Kuznetsov [30]. General solutions for equation (17) and (19) are given, respectively, by

$$u_1\left(y\right)=-{\displaystyle \frac{y^2}{2}}+c_{1}y+d_1;\eta <y<1$$

(25)

$${u_3\left(y\right)=c}_3\exp \left({\lambda }_{3}y\right)+d_3\exp \left(-{\lambda }_{3}y\right)+{\displaystyle \frac{1}{M_3{\lambda }_3^2}};0<y<\xi$$

(26)

Upon using the transformation $\tilde{y}=(\eta -y){\tau }_n$ and writing $u_2\left(y\right)\equiv U_2\left(\tilde{y}\right)$, equation (18) can be written as:

$${\displaystyle \frac{d^2{U}_2\left(\tilde{y}\right)}{d{\tilde{y}}^2}}-{\left(\tilde{y}\right)}^n{U}_2\left(\tilde{y}\right)+{\displaystyle \frac{1}{M_2{\left({\tau }_n\right)}^2}}=0$$

(27)

Equation (27) is recognized as the generalized inhomogeneous Airy’s differential equation, Swanson and Headley [42]. A fundamental set of linearly independent solutions for the homogeneous part of (27) are the generalized Airy’s functions $A_n\left(\tilde{y}\right)$ and $B_n\left(\tilde{y}\right)$, defined respectively by, [42]:

$$A_n\left(\tilde{y}\right)={\displaystyle \frac{2p}{\pi }}\sin \left(p\pi \right){\left(\tilde{y}\right)}^{1/2}{K}_p\left(\zeta \right)$$

(28)

$$B_n\left(\tilde{y}\right)={\left(p\tilde{y}\right)}^{1/2}[I_p\left(\zeta \right)+I_{-p}\left(\zeta \right)]$$

(29)

The Wronskian of $A_n\left(\tilde{y}\right)$ and $B_n\left(\tilde{y}\right)$ is given by:

$$W\left(A_n\left(\tilde{y}\right),B_n\left(\tilde{y}\right)\right)={\displaystyle \frac{2}{\pi }}{p}^{{\displaystyle \frac{1}{2}}}\sin \left(p\pi \right)$$

(30)

and the terms $I_{p}and{K}_p$ are the modified Bessel functions defined as:

$$I_p\left(\zeta \right)={\left(i\right)}^{-p}{K}_p\left(j\zeta \right)=\sum _{m=1}^{\infty }{\displaystyle \frac{1}{m!\Gamma (m+p+1)}}({\displaystyle \frac{\zeta }{2}}{)}^{2m+p}$$

(31)

$$K_p\left(\zeta \right)={\displaystyle \frac{\pi }{2}}{\displaystyle \frac{\left(I_{-p}\right(\zeta )-I_p(\zeta \left)\right)}{sin\left(p\pi \right)}}$$

(32)

with $i=\sqrt{-1}$ , $p={\displaystyle \frac{1}{n+2}}$, $\zeta =2p{\tilde{y}}^{{\displaystyle \frac{1}{2p}}}$, and $\Gamma (.)$ is the gamma function. Solution to the homogeneous part of (27) thus takes the form

$$U_{2h}\left(\tilde{y}\right)=c_2{A}_n\left(\tilde{y}\right)+d_2{B}_n\left(\tilde{y}\right)$$

(33)

A particular solution to equation (27) is expressed in the form

U2𝑝y = σnNn(y)

(34)

wherein

$${\sigma }_n={\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}$$

(35)

$$N_n\left(\tilde{y}\right)=A_n\left(\tilde{y}\right){\int }_0^{\tilde{y}}{B}_n\left(t\right)dt-B_n\left(\tilde{y}\right){\int }_0^{\tilde{y}}{A}_n\left(t\right)dt$$

(36)

The function $N_n\left(\tilde{y}\right)$, defined in (36), is known as the Generalized Nield-Kuznetsov Function of the First Kind, [34]. When $n=1$, equation (27) reduces to Airy’s equation and (36) reduces to the Standard Nield-Kuznetsov Function of the First Kind, $N_i\left(\tilde{y}\right)$, defined by Nield and Kuznetsov in [30].

General solution to equation (27) thus takes the form

$$U_2\left(\tilde{y}\right)=c_2{A}_n\left(\tilde{y}\right)+d_2{B}_n\left(\tilde{y}\right)+{\sigma }_n{N}_n\left(\tilde{y}\right)$$

(37)

Upon using $\tilde{y}=(\eta -y){\tau }_n$ and $u_2\left(y\right)\equiv U_2\left(\tilde{y}\right)$, general solution (37) is expressed as:

$$u_2\left(y\right)=c_2{A}_n\left((\eta -y){\tau }_n\right)+d_2{B}_n\left((\eta -y){\tau }_n\right)+{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}{N}_n\left(\left(\eta -y\right){\tau }_n\right)$$

(38)

Derivatives of the functions $A_n\left(\tilde{y}\right)$, $B_n\left(\tilde{y}\right)$ and $N_n\left(\tilde{y}\right)$ are given by:

$${A\prime }_n\left(\tilde{y}\right)=-{\displaystyle \frac{2p}{\pi }}\sin \left(p\pi \right){\left(\tilde{y}\right)}^{{\displaystyle \frac{n+1}{2}}}{K}_{p-1}(\zeta )$$

(39)

$${B\prime }_n\left(\tilde{y}\right)=p^{1/2}{\tilde{y}}^{{\displaystyle \frac{n+1}{2}}}[I_{p-1}\left(\zeta \right)+I_{1-p}\left(\zeta \right)]$$

(40)

$${N\prime }_n\left(\tilde{y}\right)={A\prime }_n\left(\tilde{y}\right){\int }_0^{\tilde{y}}{B}_n\left(t\right)dt-{B\prime }_n\left(\tilde{y}\right){\int }_0^{\tilde{y}}{A}_n\left(t\right)dt$$

(41)

Upon using boundary and interfacial conditions, (22)-(24), the constants in equations (25), (26) and (38), satisfy the following system of linear equations, written in the form $\mathbf{A}\mathbf{X}=\mathbf{C}$, where:

$$A=\left(\begin{array}{cccccc}1& 1& 0& 0& 0& 0\\ \eta & 1& -A_n\left(0\right)& -B_n\left(0\right)& 0& 0\\ 1& 0& {\displaystyle \frac{{\mu }_{e2}}{\mu }}{\tau }_n{A^{\prime }}_n\left(0\right)& {\displaystyle \frac{{\mu }_{e2}}{\mu }}{\tau }_n{B^{\prime }}_n\left(0\right)& 0& 0\\ 0& 0& A_n\left({\tau }_n(\eta -\xi )\right)& B_n\left({\tau }_n(\eta -\xi )\right)& -\exp \left({\lambda }_3\xi \right)& -\exp \left(-{\lambda }_3\xi \right)\\ 0& 0& {{\tau }_n{A}^{\prime }}_n\left({\tau }_n\left(\eta -\xi \right)\right)& {\tau }_n{B^{\prime }}_n\left({\tau }_n\left(\eta -\xi \right)\right)& {\displaystyle \frac{{\mu }_{e3}}{{\mu }_{e2}}}{\lambda }_3\exp \left({\lambda }_3\xi \right)& {-\lambda }_3{\displaystyle \frac{{\mu }_{e3}}{{\mu }_{e2}}}\exp \left(-{\lambda }_3\xi \right)\\ 0& 0& 0& 0& 1& 1\end{array}\right)$$

(42)

$$X=\left[\begin{array}{c}\begin{array}{c}{c}_1\\ d_1\\ c_2\end{array}\\ d_2\\ c_3\\ d_3\end{array}\right];C=\left[\begin{array}{c}\begin{array}{c}{\displaystyle \frac{1}{2}}\\ {\displaystyle \frac{{\eta }^2}{2}}+{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}{N}_n\left(0\right)\\ \eta -{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\tau }_n}}{\displaystyle \frac{{\mu }_{e2}}{\mu }}{N\prime }_n\left(0\right)\}\end{array}\\ {\displaystyle \frac{1}{M_3{\lambda }_3^2}}-{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}{N}_n\left({\tau }_n\left(\eta -\xi \right)\right).\\ -{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\tau }_n}}{N^{\prime }}_n\left({\tau }_n\left(\eta -\xi \right)\right)\\ -{\displaystyle \frac{1}{M_3{\lambda }_3^2}}\end{array}\right]$$

(43)

In solving the above linear system, we take ${\displaystyle \frac{{\mu }_{e2}}{\mu }}=$ ${\displaystyle \frac{{\mu }_{e3}}{{\mu }_{e2}}}=1$ and make use of the following values of the generalized functions $A_n$, $B_n$ and $N_n$, and their first derivatives at zero, namely:

$$\begin{gathered} N_n\left(0\right)={N^{\prime }}_n\left(0\right)=0;A_n\left(0\right)={\displaystyle \frac{{\left(p\right)}^{1-p}}{\mathsf{\Gamma }(1-p)}};B_n\left(0\right)={\displaystyle \frac{{\left(p\right)}^{1/2-p}}{\mathsf{\Gamma }(1-p)}};{A^{\prime }}_n\left(0\right)= \\ -{\displaystyle \frac{{\left(p\right)}^p}{\mathsf{\Gamma }\left(p\right)}};{B^{\prime }}_n\left(0\right)=-{\displaystyle \frac{{\left(p\right)}^{p-\frac{1}{2}}}{\mathsf{\Gamma }\left(p\right)}}. \end{gathered}$$

Following Swanson and Headley, [42], the generalized Airy’s functions are evaluated using the following relationships:

$${\rho }_n={\displaystyle \frac{({p)}^{1-p}}{\mathsf{\Gamma }(1-p)}}\mathrm{and}{\phi }_n={\displaystyle \frac{{\left(p\right)}^p}{\mathsf{\Gamma }\left(p\right)}}$$

(44)

$$g_{n1}\left(y\right)=1+{\sum }_{k=1}^{\infty }{p}^{2k}\prod _{j=1}^k{\displaystyle \frac{x^{\left(n+2\right)k}}{j\left(j-p\right)}}$$

(45)

$$g_{n2}\left(y\right)=y\left[1+{\sum }_{k=1}^{\infty }{p}^{2k}\prod _{j=1}^k{\displaystyle \frac{x^{\left(n+2\right)k}}{j\left(j+p\right)}}\right]$$

(46)

$$A_n\left(y\right)={\rho }_n{g}_{n1}\left(y\right)-{\phi }_n{g}_{n2}\left(y\right)$$

(47)

$$B_n\left(y\right)={{\displaystyle \frac{1}{\sqrt{p}}}[\rho }_n{g}_{n1}\left(y\right)+{\phi }_n{g}_{n2}\left(y\right)]$$

(48)

and the Generalized Nield-Kuznetsov Function of the First Kind can be evaluated using the following expression:

$$N_n\left(y\right)={\displaystyle \frac{2}{\sqrt{p}}}{\rho }_n{\phi }_n\left\{g_{n1}\left(y\right){\int }_0^y{g}_{n2}\left(t\right)dt-g_{n2}\left(y\right){\int }_0^y{g}_{n1}\left(t\right)dt\right\}$$

(49)

3. Results and Discussion

Results have been obtained for a range of values of n and the range of Da = 1; 0.1; 0.001; 0.0001; 0.00001; and 0.000001 in order to illustrate the effects of changing n and Da on the generalized functions, on the permeability function, on the arbitrary constants, and on the velocity profiles. Thick and thin transition layers are considered. In particular, we choose $\xi ={\displaystyle \frac{1}{3}}$ , $\eta ={\displaystyle \frac{2}{3}}$ for thick transition layer and $\xi =0.49$ , $\eta =0.51$ for thin transition layer.

3.1. Permeability Distributions

Dimensionless permeability distribution in the transition layer is given by equation (12), namely ${\displaystyle \frac{K_2}{H^2}}=Da{\left({\displaystyle \frac{\eta -\xi }{\eta -y}}\right)}^n$, and its dependence on the values of n, Da, and the thickness of the porous layer are illustrated in Figure 2, Figure 3 and Figure 4, for the thick and thin transition layer, respectively, Da = 1, and various values of n. These figures demonstrate how n acts as a control parameter in the permeability distribution. Figure 4 shows the effect of Da on the permeability distribution for n = 1 and thin transition layer. The graph is zoomed in to better illustrate the distribution when Da = 0.000001.

3.2. Velocity Profiles

Figure 5 illustrates the velocity profiles across the three layer flow configuration for the case of thin transition layer, n =1, and various values of Da. The transition layer thickness is shown on the graphs as a small thick line (smudge) on each graph.

In using the solution procedure given by equations (44)-(49), the smallest Darcy number we could compute results for was $Da=0.0000005.$ However, results are shown for a smallest value of $Da=0.000001$, when $n=1$. Figure 5 shows smooth transitions from one flow domain layer to another for the range of $Da$ used. As Da increases, the singularity at the lower interface starts creeping in (as anticipated by Nield and Kuznetsov [ ]). When n=2, Figure 6 illustrates a similar pattern of velocity profiles, and shows a more noticeable effect of the singularity at the lower interface. This emphasizes that the current model and solution procedure is more suitable for low values of Da.

For different values of n, and a fixed value of Da, the velocity profiles are close to each other (due to graph scale) and follow a similar pattern and shape to the profiles of Figure 6. Therefore, to better demonstrate the relative magnitudes of velocities for various values of n, velocity profiles are sketched over the span of the thin transition layer $0.49\le \eta -\xi \le 0.51$, illustrated in Figure 7. The increase in permeability in the transition layer with increasing n is accompanied with an increase in velocity, as Figure 7 demonstrates.

3.3. Velocity at the Interfaces

At the lower and upper interfaces, $y=\xi$ and $y=\eta$, respectively, velocity expressions are obtained from equations (25), (26), and (38), and take the form

$$u_1\left(\eta \right)=-{\displaystyle \frac{{\eta }^2}{2}}+c_1\eta +d_1$$

(50)

$$u_2\left(\eta \right)=c_2{A}_n\left(0\right)+d_2{B}_n\left(0\right)$$

(51)

$${u_3\left(\xi \right)=c}_3\mathit{exp}\left({\lambda }_3\xi \right)+d_3\mathit{exp}\left(-{\lambda }_3\xi \right)+{\displaystyle \frac{1}{M_3{\lambda }_3^2}}$$

(52)

$$u_2\left(\xi \right)=c_2{A}_n\left((\eta -\xi ){\tau }_n\right)+d_2{B}_n\left((\eta -\xi ){\tau }_n\right)+{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}{N}_n(\left(\eta -\xi \right){\tau }_n)$$

(53)

Velocity at the lower and upper interfaces for thick and thin transition layers, for different values of Da and different values of n, are given in Table 2 and Table 3. In Table 2, velocity at the upper interface is illustrated for the thick transition layer, only, since the behaviour is similar for the thin layer, for n=1,2,3 and various values of Da.

As n increases, for a fixed value of Da, interfacial velocity increases since permeability increases as we get closer to the upper interface. This behaviour is consistent for all values of Da shown. Interfacial velocity consistently increases with increasing Da for any fixed value of n, as expected.

In Table 3, velocity at the lower interface is illustrated for the thin transition layer, only, since the behaviour is similar for the thick layer, for n=1,2,3 and various values of Da.

As n increases, for a fixed value of Da, interfacial velocity decreases. This is due to the velocity continuity condition at the lower interface where the velocity in the homogeneous, constant permeability layer has the effect of slowing down the flow velocity at the interface.

This behaviour is consistent for all values of Da shown. Interfacial velocity consistently increases with increasing Da for any fixed value of n, since Da is defined in terms of a fixed, reference permeability.

Agreement between the values computed using equations (50) and (51) was for more than 12 significant digits, which indicates that the current formulation and computations are accurate. Similarly when using equations (52) and (53). In Table 2 and Table 3 we retained the first 8 digits for illustration.

3.4. Shear Stress at the Interfaces

At the lower and upper interfaces, $y=\xi$ and $y=\eta$, respectively, shear stress expressions are obtained from equations (25), (26), and (38), and take the form

$${\displaystyle \frac{d{u}_1}{dy}}\left(\eta \right)=c_1-\eta$$

(54)

$${\displaystyle \frac{d{u}_2}{dy}}\left(\eta \right)={-c}_2{{\tau }_{n}A\prime }_n\left(0\right)-d_2{{\tau }_n{B}^{\prime }}_n\left(0\right)$$

(55)

$${{\displaystyle \frac{d{u}_3}{dy}}\left(\xi \right)=c}_3{\lambda }_3\exp \left({\lambda }_3\xi \right)-d_3{\lambda }_3\exp \left(-{\lambda }_3\xi \right)$$

(56)

$$\begin{gathered} {\displaystyle \frac{{du}_2}{dy}}\left(\xi \right)={-c}_2{\tau }_n{A\prime }_n\left((\eta -\xi ){\tau }_n\right)-d_2{{\tau }_n{B}^{\prime }}_n\left(\left(\eta -\xi \right){\tau }_n\right)- \\ {\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\tau }_n}}{N\prime }_n\left(\left(\eta -\xi \right){\tau }_n\right). \end{gathered}$$

(57)

Shear stress values at the lower and upper interfaces for different middle layer thickness, and for n = 1, 2, 3, are given in Table 4, which furnish the following observations.

• For a given Da, the value of the shear stress at the upper interface decreases with increasing n for the thick layer. However, it increases with increasing n for thin layer.
• For a given Da, the value of the shear stress at the upper interface increases with increasing n for the thin layer, and decreases with increasing n for the thick layer.
• The thin layer experiences a higher sheer than the thick layer at the upper interface with increasing n.
• The thick layer experiences a higher sheer than the thin layer at the lower interface when n deviates from 1.

3.5. Mean Velocity Across the Layers

The dimensionless mean velocities across the lower and upper layers are defined, respectively, as

$${\overline{u}}_1={\int }_{\eta }^1{u}_{1}dy={\displaystyle \frac{1}{6}}\left({\eta }^3-1\right)-{\displaystyle \frac{c_1}{2}}\left({\eta }^2-1\right)-d_1(\eta -1)$$

(58)

$${\overline{u}}_3={\int }_0^{\xi }{u}_{3}dy={\displaystyle \frac{c_3}{{\lambda }_3}}\left(e^{{\lambda }_3\xi }-1\right)-{\displaystyle \frac{d_3}{{\lambda }_3}}\left(e^{{-\lambda }_3\xi }-1\right)+{\displaystyle \frac{\xi }{M_3{\lambda }_3^2}}$$

(59)

Dimensionless mean velocity in the transition layer is obtained from (38) and can be written as:

$${\overline{u}}_2={\int }_{\xi }^{\eta }{u}_{2}dy=-{\displaystyle \frac{1}{{\tau }_n}}{\int }_{(\eta -\xi ){\tau }_n}^0\left(c_2{A}_n\left(t\right)+d_2{B}_n\left(t\right)+{\displaystyle \frac{\pi }{2\sqrt{p}\sin \left(p\pi \right)M_2{\left({\tau }_n\right)}^2}}{N}_n\left(t\right)\right)dt$$

(60)

where $t=(\eta -y){\tau }_n$.

The dimensionless mean velocity across the configuration (that is, across the three layers) is defined as

$$\overline{u}={\overline{u}}_1+{\overline{u}}_2+{\overline{u}}_3$$

(61)

and represents a measure of the overall volume flux through the channel configuration.

The above expressions are evaluated using Maple and the values are listed in for n = 1 and n = 2 in Table 5 for different transition layer thicknesses and Da = 1, 0.1, and 0.01 so that we can compare our results with those obtained by Nield and Kuznetsov [30]. To make the comparison, we reformulated and solved the flow problem using the same configuration as that of Nield and Kuznetsov [30]. Agreement between the results is evident from Table 5.

For a given transition layer thickness and a given value of n, the total volume flux decreases with decreasing Da. This is due to flow retardation for smaller Da. For a given Da, the total volume flux increases with increasing n.

4. Conclusions

In this work, we considered the problem of flow through a free-space channel over a porous layer in the presence of a transition porous layer. This problem was treated by Nield and Kuznetsov [30] to illustrate the characteristics of the flow when the transition layer is a Brinkman layer of variable permeability described as the reciprocal of a linear polynomial. In this work we introduced a generalized form of the permeability function that involves the reciprocal of the nth degree of the a linear polynomial. This choice of permeability function resulted in the inhomogeneous, generalized Airy’s differential equation. We provided the solution to the resulting differential equation in terms of a generalized Nield-Kuznetsov function, together with a computational procedure to evaluate the arising special functions. For the case of n = 1, the current problem reduces to the one considered by Nield and Kuznetsov, [30], and our results agree with their results for total volume flux across the tri-layer configuration.

Highlights of this work may be stated as follows:

• The introduction of a generalized permeability function that provides modeling flexibility and validity for small Darcy number, and possible control over permeability amplification in the transition layer.
• Obtaining the general and particular solutions to the resulting inhomogeneous generalized Airy’s equation through the introduction of a generalized Nield-Kuznetsov function.
• Providing an evaluation procedure for the arising generalised Airy’s functions and the generalized Nield-Kuznetsov function.