In the above equations, $G=-\frac{\text{d}p}{\text{d}x}$ is the constant pressure gradient. For

$i=1,2,3$ , with $i=1$ referring to layer 1, $i=2$ referring to layer 2, and $i=3$ referring to layer 3, the quantities in Equations (4), (5), and (6): ${u}_{i}^{\mathrm{*}}\mathrm{,}{K}_{i}\mathrm{,}{\mu}_{i}$ , and ${\mu}_{ei}$ denote the velocity, permeability, base fluid viscosity, and effective viscosity, respectively, in layer $i$ .

Following Nield and Kozentsov, [1] , we introduce the following dimensionless

variables, where $Da$ refers to Darcy number: $y=\frac{{y}^{*}}{H};Da=\frac{{K}_{0}}{{H}^{2}};{u}_{i}=\frac{{\mu}_{i}}{G{H}^{2}}{u}_{i}^{*}$ for $i=1,2,3$ ; and ${M}_{i}=\frac{{\mu}_{ei}}{{\mu}_{i}}$ for $i=2,3$ .

The permeability distributions, ${K}_{1},{K}_{2}$ and ${K}_{3}$ can be written in the following dimensionless form in terms of Darcy number, Da:

${K}_{1}^{*}=\frac{{K}_{1}}{{H}^{2}}\to \infty $ (7)

${K}_{2}^{*}=\frac{{K}_{2}}{{H}^{2}}=Da{\left(\frac{\xi -\eta}{y-\eta}\right)}^{n}$ (8)

${K}_{3}^{*}=\frac{{K}_{3}}{{H}^{2}}=Da$ . (9)

Equations (4), (5), and (6) take the following forms, respectively:

$\frac{{\text{d}}^{2}{u}_{1}}{\text{d}{y}^{2}}+1=0;\text{\hspace{1em}}0<y<\eta .$ (10)

${M}_{2}\frac{{\text{d}}^{2}{u}_{2}}{\text{d}{y}^{2}}-\frac{{H}^{2}}{{K}_{2}}{u}_{2}+1=0;\text{\hspace{1em}}\eta <y<\xi .$ (11)

${M}_{3}\frac{{\text{d}}^{2}{u}_{3}}{\text{d}{y}^{2}}-\frac{{H}^{2}}{{K}_{3}}{u}_{3}+1=0;\text{\hspace{1em}}\xi <y<1$ (12)

and, upon substituting the permeability distributions, (7)-(9), in equations (10)-(12), we get

$\frac{{\text{d}}^{2}{u}_{1}}{\text{d}{y}^{2}}+1=0;\text{\hspace{1em}}0<y<\eta .$ (13)

$\frac{{\text{d}}^{2}{u}_{2}}{\text{d}{y}^{2}}-\frac{1}{{M}_{2}Da{\left(\xi -\eta \right)}^{n}}{\left(y-\eta \right)}^{n}{u}_{2}+\frac{1}{{M}_{2}}=0;\text{\hspace{1em}}\eta <y<\xi .$ (14)

$\frac{{\text{d}}^{2}{u}_{3}}{\text{d}{y}^{2}}-\frac{1}{Da{M}_{3}}{u}_{3}+\frac{1}{{M}_{3}}=0;\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\xi <y<1.$ (15)

Equation (13), (14) and (15) must be solved subject to the following boundary and matching conditions

${u}_{1}=0\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}y=0$ (16a)

${u}_{1}={u}_{2}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}y=\eta $ (16b)

$\frac{\text{d}{u}_{1}}{\text{d}y}={\vartheta}_{1}\frac{\text{d}{u}_{2}}{\text{d}y}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}\text{\hspace{0.05em}}y=\eta .$ (16c)

${u}_{2}={u}_{3}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}y=\xi .$ (16d)

${\vartheta}_{2}\frac{\text{d}{u}_{2}}{\text{d}y}=\frac{\text{d}{u}_{3}}{\text{d}y}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}\text{\hspace{0.05em}}y=\xi .$ (16e)

${u}_{3}=0\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{at}\text{\hspace{0.17em}}y=1$ , (16f)

where ${\vartheta}_{1}=\frac{{\mu}_{2}}{{\mu}_{1}}$ and ${\vartheta}_{2}=\frac{{\mu}_{2}}{{\mu}_{3}}$ are the viscosity ratios.

Now, letting

${\tau}_{n}=\frac{1}{\sqrt[n+2]{{M}_{2}Da{\left(\xi -\eta \right)}^{n}}}$ (17)

${\lambda}_{3}=\frac{1}{\sqrt{{M}_{3}Da}}$ (18)

Equation (14) and Equation (15) become, respectively:

$\frac{{\text{d}}^{2}{u}_{2}}{\text{d}{y}^{2}}-{\left({\tau}_{n}\right)}^{n+2}{\left(y-\eta \right)}^{n}+\frac{1}{{M}_{2}}=0$ (19)

$\frac{{\text{d}}^{2}{u}_{3}}{\text{d}{y}^{2}}-{\lambda}_{3}^{2}{u}_{3}+\frac{1}{{M}_{3}}=0$ . (20)

General solutions for Equation (13) and Equation (20) are given, respectively, by

${u}_{1}\left(y\right)={c}_{1}y+{d}_{1}-\frac{{y}^{2}}{2}$ (21)

${u}_{3}\left(y\right)={c}_{3}\mathrm{exp}\left({\lambda}_{3}y\right)+{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}y\right)+\frac{1}{{M}_{3}{\lambda}_{3}^{2}}$ . (22)

In order to solve Equation (19), we first use the following transformation $\stackrel{\u02dc}{y}={\tau}_{n}\left(y-\eta \right)$ and write ${u}_{2}\left(y\right)\equiv {U}_{2}\left(\stackrel{\u02dc}{y}\right)$ . Equation (19) then becomes:

$\frac{{\text{d}}^{2}{U}_{2}\left(\stackrel{\u02dc}{y}\right)}{\text{d}{\stackrel{\u02dc}{y}}^{2}}-{\left(\stackrel{\u02dc}{y}\right)}^{n}{U}_{2}\left(\stackrel{\u02dc}{y}\right)+\frac{1}{{M}_{2}{\left({\tau}_{n}\right)}^{2}}=0$ (23)

Equation (23) is the generalized inhomogeneous Airy’s differential equation. A fundamental pair of linearly independent solutions for the homogeneous part are the generalized Airy’s functions ${A}_{n}\left(\stackrel{\u02dc}{y}\right)$ and ${B}_{n}\left(\stackrel{\u02dc}{y}\right)$ , (cf. Swanson and Headley [13] ), defined respectively by

${A}_{n}\left(\stackrel{\u02dc}{y}\right)=\frac{2p}{\text{\pi}}\mathrm{sin}\left(p\text{\pi}\right){\left(\stackrel{\u02dc}{y}\right)}^{1/2}{K}_{p}\left(\zeta \right)$ (24)

${B}_{n}\left(\stackrel{\u02dc}{y}\right)={\left(p\stackrel{\u02dc}{y}\right)}^{\frac{1}{2}}\left({I}_{-p}\left(\zeta \right)+{I}_{p}\left(\zeta \right)\right)$ (25)

The terms ${I}_{p}$ and ${K}_{p}$ are the modified Bessel functions defined as:

${I}_{p}\left(\zeta \right)={i}^{-p}{J}_{p}\left(i\zeta \right)={\displaystyle {\sum}_{m=1}^{\infty}\frac{1}{m!\Gamma \left(m+p+1\right)}{\left(\frac{\zeta}{2}\right)}^{2m+p}}$ (26)

${K}_{p}(\zeta )=\frac{\text{\pi}}{2}\frac{\left({I}_{-p}\left(\zeta \right)-{I}_{p}\left(\zeta \right)\right)}{\mathrm{sin}\left(p\text{\pi}\right)}$ (27)

with $p=\frac{1}{n+2}$ , $\zeta =2p{Y}^{\frac{1}{2p}}$ , and $\Gamma $ is the gamma function.

Solution to the homogeneous part of (23) is thus given by

${U}_{2h}\left(\stackrel{\u02dc}{y}\right)={c}_{2}{A}_{n}\left(\stackrel{\u02dc}{y}\right)+{d}_{2}{B}_{n}\left(\stackrel{\u02dc}{y}\right)$ . (28)

We find it convenient to introduce the following integral function:

${Z}_{n}\left(\stackrel{\u02dc}{y}\right)={A}_{n}\left(\stackrel{\u02dc}{y}\right){\displaystyle {\int}_{0}^{\stackrel{\u02dc}{y}}{B}_{n}\left(t\right)\text{d}t}-{B}_{n}\left(\stackrel{\u02dc}{y}\right){\displaystyle {\int}_{0}^{\stackrel{\u02dc}{y}}{A}_{n}\left(t\right)\text{d}t}.$ (29)

The function ${Z}_{n}\left(\stackrel{\u02dc}{y}\right)$ reduces to the Nield-Kuznetsov function ${N}_{i}\left(\stackrel{\u02dc}{y}\right)$ when n = 1.

The Wronskian of ${A}_{n}\left(Y\right)$ and ${B}_{n}\left(Y\right)$ is given by:

$W\left({A}_{n}\left(\stackrel{\u02dc}{y}\right),{B}_{n}\left(\stackrel{\u02dc}{y}\right)\right)=\frac{2}{\text{\pi}}{p}^{\frac{1}{2}}\mathrm{sin}\left(p\text{\pi}\right)$ (30)

and general solution of (23) is expressed, using variation of parameters, as:

${U}_{2}\left(\stackrel{\u02dc}{y}\right)\mathrm{=}{c}_{2}{A}_{n}\left(\stackrel{\u02dc}{y}\right)+{d}_{2}{B}_{n}\left(\stackrel{\u02dc}{y}\right)+{\sigma}_{n}{Z}_{n}\left(\stackrel{\u02dc}{y}\right)$ (31)

where $p=\frac{1}{n+2}$ and ${\sigma}_{n}=\frac{\text{\pi}}{2\sqrt{p}sin\left(p\text{\pi}\right){M}_{2}{\left({\tau}_{n}\right)}^{2}}$ .

Upon substituting $\stackrel{\u02dc}{y}={\tau}_{n}\left(y-\eta \right)$ , ${u}_{2}\left(y\right)\equiv {U}_{2}\left(\stackrel{\u02dc}{y}\right)$ , and ${\sigma}_{n}=\frac{\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\left({\tau}_{n}\right)}^{2}}$ , in (28) we obtain the following general solution to Equation (14):

$\begin{array}{c}{u}_{2}\left(y\right)={c}_{2}{A}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)+{d}_{2}{B}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}sin\left(p\text{\pi}\right){M}_{2}{\left({\tau}_{n}\right)}^{2}}{Z}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)\end{array}$ (32)

Derivatives of the functions ${A}_{n}\left(\stackrel{\u02dc}{y}\right)$ , ${B}_{n}\left(\stackrel{\u02dc}{y}\right)$ and ${Z}_{n}\left(\stackrel{\u02dc}{y}\right)$ are given by:

${{A}^{\prime}}_{n}\left(\stackrel{\u02dc}{y}\right)=-\frac{2p}{\text{\pi}}\left[\mathrm{sin}py\right]{\stackrel{\u02dc}{y}}^{\frac{n+1}{2}}{K}_{p-1}\left(\zeta \right)$ (33)

${{B}^{\prime}}_{n}\left(\stackrel{\u02dc}{y}\right)={p}^{\frac{1}{2}}{\stackrel{\u02dc}{y}}^{\frac{n+1}{2}}\left[{I}_{1-p}\left(\zeta \right)+{I}_{p-1}\left(\zeta \right)\right]$ (34)

${{Z}^{\prime}}_{n}\left(\stackrel{\u02dc}{y}\right)={{A}^{\prime}}_{n}\left(\stackrel{\u02dc}{y}\right){\displaystyle {\int}_{0}^{\stackrel{\u02dc}{y}}{B}_{n}\left(t\right)\text{d}t}-{{B}^{\prime}}_{n}\left(\stackrel{\u02dc}{y}\right){\displaystyle {\int}_{0}^{\stackrel{\u02dc}{y}}{A}_{n}\left(t\right)\text{d}t}.$ (35)

Shear stress expressions across the layers are obtained from Equations (21), (22) and (32), and take the following form:

$\frac{\text{d}{u}_{1}}{\text{d}y}={c}_{1}-y$ (36)

$\begin{array}{c}\frac{\text{d}{u}_{2}}{\text{d}y}={c}_{2}{\tau}_{n}{{A}^{\prime}}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)+{d}_{2}{\tau}_{n}{{B}^{\prime}}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\tau}_{n}}{{Z}^{\prime}}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)\end{array}$ (37)

$\frac{\text{d}{u}_{3}}{\text{d}y}={c}_{3}{\lambda}_{3}\mathrm{exp}\left({\lambda}_{3}y\right)-{d}_{3}{\lambda}_{3}\mathrm{exp}\left(-{\lambda}_{3}y\right)$ . (38)

Upon using boundary and interfacial conditions, (16a)-(16f), in Equations (21), (22) and (32), we obtain the following system of linear equations that is to be solved for the arbitrary constants ${c}_{1},{d}_{1},{c}_{2},{d}_{2},{c}_{3},{d}_{3}$ :

${d}_{1}=0$ (39)

${c}_{1}\eta +{d}_{1}-{c}_{2}{A}_{n}\left(0\right)-{d}_{2}{B}_{n}\left(0\right)=\frac{{\eta}^{2}}{2}$ (40)

${c}_{1}-{c}_{2}{\vartheta}_{1}{\tau}_{n}{{A}^{\prime}}_{n}\left(0\right)-{d}_{2}{\vartheta}_{1}{\tau}_{n}{{B}^{\prime}}_{n}\left(0\right)=\eta $ (41)

$\begin{array}{l}{c}_{2}{A}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)+{d}_{2}{B}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)-{c}_{3}\mathrm{exp}\left({\lambda}_{3}\xi \right)-{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}\xi \right)\\ =\frac{1}{{M}_{3}{\lambda}_{3}^{2}}-\frac{\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\left({\tau}_{n}\right)}^{2}}{Z}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\end{array}$ (42)

$\begin{array}{l}{c}_{2}{\vartheta}_{2}{\tau}_{n}{{A}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)+{d}_{2}{\vartheta}_{2}{\tau}_{n}{{B}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)-{\lambda}_{3}{c}_{3}\mathrm{exp}\left({\lambda}_{3}\xi \right)+{\lambda}_{3}{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}\xi \right)\\ =-\frac{{\vartheta}_{2}\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\tau}_{n}}{{Z}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\end{array}$ (43)

${c}_{3}\mathrm{exp}\left({\lambda}_{3}\right)+{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}\right)=\frac{-1}{{M}_{3}{\lambda}_{3}^{2}}$ (44)

Linear Equations (39)-(44) are cast in the following matrix-vector form

$Mx=c$ , (45)

where:

$M=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ \eta & 1& -{A}_{n}\left(0\right)& -{B}_{n}\left(0\right)& 0& 0\\ 1& 0& -{\vartheta}_{1}{\tau}_{n}{{A}^{\prime}}_{n}\left(0\right)& -{\vartheta}_{1}{\tau}_{n}{B}_{n}\left(0\right)& 0& 0\\ 0& 0& {{A}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)& {B}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)& -\mathrm{exp}\left({\lambda}_{3}\xi \right)& -\mathrm{exp}\left(-{\lambda}_{3}\xi \right)\\ 0& 0& {\vartheta}_{2}{\tau}_{n}{{A}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)& {\vartheta}_{2}{\tau}_{n}{{B}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)& -{\lambda}_{3}\mathrm{exp}\left({\lambda}_{3}\xi \right)& {\lambda}_{3}\mathrm{exp}\left(-{\lambda}_{3}\xi \right)\\ 0& 0& 0& 0& \mathrm{exp}\left({\lambda}_{3}\right)& \mathrm{exp}\left(-{\lambda}_{3}\right)\end{array}\right)$ (46)

$x=\left[\begin{array}{c}{c}_{1}\\ {d}_{1}\\ {c}_{2}\\ {d}_{2}\\ {c}_{3}\\ {d}_{3}\end{array}\right]$ (47)

$c=\left(\begin{array}{c}0\\ \frac{{\eta}^{2}}{2}\\ \eta \\ \frac{1}{{M}_{3}{\lambda}_{3}^{2}}-\frac{\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{({\tau}_{n})}^{2}}{Z}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\\ -\frac{{\vartheta}_{2}\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\tau}_{n}}{{Z}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\\ \frac{-1}{{M}_{3}{\lambda}_{3}^{2}}\end{array}\right)$ . (48)

In solving the above linear system, we make use of the following values of the generalized functions ${A}_{n}$ , ${B}_{n}$ and ${Z}_{n}$ , and their first derivatives at zero (cf. [13] ), wherein $\Gamma $ is the Gamma function:

${A}_{n}\left(0\right)=\frac{{p}^{1-p}}{\Gamma \left(1-p\right)}$ (49)

${B}_{n}\left(0\right)=\frac{{p}^{\frac{1}{2}-p}}{\Gamma \left(1-p\right)}$ (50)

${{A}^{\prime}}_{n}\left(0\right)=-\frac{{p}^{p}}{\Gamma \left(p\right)}$ (51)

${{B}^{\prime}}_{n}\left(0\right)=\frac{{p}^{p-\frac{1}{2}}}{\Gamma \left(p\right)}$ (52)

${Z}_{n}\left(0\right)=0$ (53)

${{Z}^{\prime}}_{n}\left(0\right)=0$ . (54)

3. Results and Analysis

Results have been obtained for a range of values of n and the range of Da = 1; 0.1; 0.001; 0.0001; and 0.00001 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 $\eta =1/3,\xi =2/3$ and $\eta =1/4,\xi =3/4$ for thick transition layer, and $\eta =0.49,\xi =0.51$ for thin layer.

Values of $p=\frac{1}{n+2}$ , ${A}_{n}$ , ${{A}^{\prime}}_{n}$ , ${B}_{n}$ , ${{B}^{\prime}}_{n}$ , ${Z}_{n}$ , ${{Z}^{\prime}}_{n}$ , ${\sigma}_{n}$ and ${\tau}_{n}$ at 0 and at

other specified points for different values of n and Da, have been computed using Maple. These values are important in the calculation of the arbitrary constants in the matrix-vector Equations (45), and in the determination and plotting of permeability functions and velocity profiles. Computational results are most accurate when Da = 1; 0.1; and 0.001.

4. Dimensionless Permeability Distributions

For the sake of graphs we find it more convenient to plot the reciprocal of the dimensionless permeability functions (7)-(9), namely

$\frac{1}{{K}_{1}^{*}}=0\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.05em}}0<y<\eta .$ (55)

$\frac{1}{{K}_{2}^{*}\left(y\right)}=\frac{1}{Da}{\left(\frac{y-\eta}{\xi -\eta}\right)}^{n};\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\eta <y<\xi .$ (56)

$\frac{1}{{K}_{1}^{*}}=\frac{1}{Da}\text{\hspace{1em}}\text{\hspace{0.05em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\xi <y<1.$ (57)

Dependence of the permeability distributions on the value of n is illustrated in Figure 2, which shows the reciprocal of the dimensionless permeability across

the three layers for n = 1, 2, 3 and 5, for $\eta =\frac{1}{3},\xi =\frac{2}{3}$ , and Da = 1. Figure 2

Figure 2. Dimensionless permeability distributions for $\eta =\frac{1}{3},\xi =\frac{2}{3}$ , Da = 1, and different n.

demonstrates the increase of the reciprocal permeability of the transition layer (or decrease of permeability) with increasing n.

5. Velocity at the Interfaces

At the lower and upper interfaces, $y=\eta $ and $y=\xi $ , respectively, between layers, velocity expressions are obtained from Equations (21), (22), and (32), and take the form

${u}_{1}\left(\eta \right)={c}_{1}\eta +{d}_{1}-\frac{{\eta}^{2}}{2}$ (58)

${u}_{2}\left(\eta \right)={c}_{2}{A}_{n}\left(0\right)+{d}_{2}{B}_{n}\left(0\right)$ (59)

$\begin{array}{c}{u}_{2}\left(\xi \right)={c}_{2}{A}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)+{d}_{2}{B}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}sin\left(p\text{\pi}\right){M}_{2}{({\tau}_{n})}^{2}}{Z}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\end{array}$ (60)

${u}_{3}\left(\xi \right)={c}_{3}\mathrm{exp}\left({\lambda}_{3}\xi \right)+{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}\xi \right)+\frac{1}{{M}_{3}{\lambda}_{3}^{2}}$ . (61)

Velocity at the lower and upper interfaces for different middle layer thickness, different values of Da, and for n = 1 and n = 2 are given in Table 1 and Table 2, which furnish the following observations:

Table 1. Velocity for n = 1 at the upper and lower interface.

Table 2. Velocity for n = 2 at the upper and lower interface.

1) Computations render reasonable results for Da as low as 0.001 when n = 1. Inaccuracies start creeping when Da < 0.001. For n = 2, results are reasonable for Da as low as 0.0001 and inaccuracies creep when Da < 0.0001. This behavior may be attributed to both round-off errors for small Da and inaccuracies in evaluation of Airy’s functions for small Da. This behavior persists for thick layers, and is less noticeable for thin middle layer calculations.

2) Computations of velocity at the lower interface using expressions (58) and (59) agree up to within seven significant digits. The same is true for computations of velocity at the upper interface using expressions (60) and (61). This is indicative of appropriate matching of the velocity at the interfaces, used in this work.

3) For a given Da, velocity at the lower interface decreases with increasing middle layer thickness. Similarly, velocity at the upper interface decreases with increasing middle layer thickness. This behavior persists for both n = 1 and n = 2.

4) For a given middle layer thickness, velocity at each of the lower and upper interfaces increases with increasing Da. This is expected, as Da is defined as a dimensionless reference permeability and accompanied with increasing permeability is a velocity increase.

5) The effect of increasing n on the velocity at the interfaces for a given thickness and Da is that the velocity at each interface increases with increasing n. This is true for both thin and thick transition layers, and for the range of Da used.

6. Shear Stress at the Interfaces

At the lower and upper interfaces, $y=\eta $ and $y=\xi $ , respectively, between layers, shear stress expressions are obtained from Equations (36), (37), and (38), and take the form

$\frac{\text{d}{u}_{1}\left(\eta \right)}{\text{d}y}={c}_{1}-\eta $ (62)

$\frac{\text{d}{u}_{2}\left(\eta \right)}{\text{d}y}={c}_{2}{\tau}_{n}{{A}^{\prime}}_{n}\left(0\right)+{d}_{2}{\tau}_{n}{{B}^{\prime}}_{n}\left(0\right)$ (63)

$\begin{array}{c}\frac{\text{d}{u}_{2}\left(\xi \right)}{\text{d}y}={c}_{2}{\tau}_{n}{{A}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)+{d}_{2}{\tau}_{n}{{B}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}\mathrm{sin}\left(p\text{\pi}\right){M}_{2}{\tau}_{n}}{{Z}^{\prime}}_{n}\left({\tau}_{n}\left(\xi -\eta \right)\right)\end{array}$ (64)

$\frac{\text{d}{u}_{3}\left(\xi \right)}{\text{d}y}={c}_{3}{\lambda}_{3}\mathrm{exp}\left({\lambda}_{3}\xi \right)-{d}_{3}{\lambda}_{3}\mathrm{exp}\left(-{\lambda}_{3}\xi \right)$ . (65)

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

1) Computations render reasonable results for Da as low as 0.001 when n = 1 and n = 2. Inaccuracies start creeping when Da < 0.001. This behavior may be attributed to the inaccuracies reported earlier when computing velocity at the interfaces using Maple’s built-in functions. When dealing with a thin transition layer, results are accurate for as low as Da = 0.00001.

2) Computations of shear stress at the lower interface using expressions (62) and (63) agree up to within a minimum of five significant digits. The same is true for computations of velocity at the upper interface using expressions (64) and (65).

3) For a given Da, the absolute value of the shear stress at the lower interface increases with increasing transition layer thickness. Similarly, at the upper interface. This behavior persists for both n = 1 and n = 2.

Table 3. Shear stress at the upper and lower interfaces for n = 1.

Table 4. Shear stress at the upper and lower interfaces for n = 2.

4) For a given transition layer thickness, absolute value of shear stress at each of the lower and upper interfaces increases with increasing Da.

5) The effect of increasing n on the absolute value of shear stress at the interfaces for a given thickness and Da is that at each interface, this absolute value increases with increasing n. This is true for both thin and thick transition layers, and for the range of Da used.

7. Friction Factor

A quantity of interest is the negative of the shear stress term at the interface

between the channel and the transition layer, namely ${c}_{f}=-\frac{\text{d}{u}_{1}}{\text{d}y}$ at $y=\eta $ .

This has been analyzed and defined by Nield and Kuznetsov, [1] , as a friction factor representing the dimensionless frictional stress in the fluid at the interface. From Equation (62) we can see that

${c}_{f}=-\frac{\text{d}{u}_{1}\left(\eta \right)}{\text{d}y}=\eta -{c}_{1}$ . (66)

Values of $-{c}_{f}$ for different Da, layer thickness and n = 1 are listed in Table 5 and compared with the values obtained by Nield and Koznetsov, [1] , for the same problem. Agreement between the current results and Nield and Koznetsov results, up to the three significant digits they report is clear from Table 5. Values of $-{c}_{f}$ for different Da, layer thickness and n = 2 are listed in Table 6.

8. Mean Velocity across the Layers

The dimensionless mean velocities across the layers are defined as

${\stackrel{\xaf}{u}}_{1}={\displaystyle \underset{0}{\overset{\eta}{\int}}{u}_{1}\text{d}y}=\frac{1}{2}{c}_{1}{\eta}^{2}+{d}_{1}\eta -\frac{1}{6}{\eta}^{3}$ (67)

$\begin{array}{c}{\stackrel{\xaf}{u}}_{2}={\displaystyle \underset{\eta}{\overset{\xi}{\int}}{u}_{2}}\text{d}y={\displaystyle \underset{\eta}{\overset{\xi}{\int}}{c}_{2}{A}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)+{d}_{2}{B}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}sin\left(p\text{\pi}\right){M}_{2}{({\tau}_{n})}^{2}}{Z}_{n}\left({\tau}_{n}\left(y-\eta \right)\right)\text{d}y\end{array}$ (68)

Table 5. Values of ${c}_{f}$ for different Da, layer thickness and n = 1.

*Nield and koznetsov results.

Table 6. Values of ${c}_{f}$ for different Da, layer thickness and n = 2.

By letting $t={\tau}_{n}\left(y-\eta \right)$ in (7.61) we get

$\begin{array}{c}{\stackrel{\xaf}{u}}_{2}=\frac{{c}_{2}}{{\tau}_{n}}{\displaystyle \underset{0}{\overset{{\tau}_{n}\left(y-\eta \right)}{\int}}{A}_{n}\left(t\right)\text{d}t}+\frac{{d}_{2}}{{\tau}_{n}}{\displaystyle \underset{0}{\overset{{\tau}_{n}\left(y-\eta \right)}{\int}}{B}_{n}\left(t\right)\text{d}t}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\text{\pi}}{2\sqrt{p}sin\left(p\text{\pi}\right){M}_{2}{\left({\tau}_{n}\right)}^{3}}{\displaystyle \underset{0}{\overset{{\tau}_{n}\left(y-\eta \right)}{\int}}{Z}_{n}\left(t\right)\text{d}t}\end{array}$ (69)

$\begin{array}{c}{\stackrel{\xaf}{u}}_{3}={\displaystyle \underset{\xi}{\overset{1}{\int}}{u}_{3}\text{d}y}={\displaystyle \underset{\xi}{\overset{1}{\int}}\left[{c}_{3}\mathrm{exp}\left({\lambda}_{3}y\right)+{d}_{3}\mathrm{exp}\left(-{\lambda}_{3}y\right)+\frac{1}{{M}_{2}{\lambda}_{3}^{2}}\right]\text{d}y}\\ =\frac{1}{{\lambda}_{3}}\left[{c}_{3}\left\{\mathrm{exp}\left({\lambda}_{3}\right)-\mathrm{exp}\left({\lambda}_{3}\xi \right)\right\}-{d}_{3}\left\{\mathrm{exp}\left(-{\lambda}_{3}\right)-\mathrm{exp}\left(-{\lambda}_{3}\xi \right)\right\}+\frac{1-\xi}{{M}_{2}{\lambda}_{3}}\right]\end{array}$ (70)

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

$\stackrel{\xaf}{u}={\stackrel{\xaf}{u}}_{1}+{\stackrel{\xaf}{u}}_{2}+{\stackrel{\xaf}{u}}_{3}={\displaystyle \underset{0}{\overset{\eta}{\int}}{u}_{1}\text{d}y}+{\displaystyle \underset{\eta}{\overset{\xi}{\int}}{u}_{2}\text{d}y}+{\displaystyle \underset{\xi}{\overset{1}{\int}}{u}_{3}}\text{d}y$ (71)

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 7 and Table 8, respectively, for different transition layer thicknesses and different Da. Comparison is made for the case of n = 1 with the Nield and Koznetsov results [1] . Agreement between the results is evident from Table 7. It is also evident from Table 7 and Table 8 that increasing n results in a higher mean velocity across the channel for thick layers but the increase is minimal for thin layers, for a given Da.

Table 7. Mean velocity across the channel for n = 1.

Table 8. Mean velocity across the channel for n = 2.

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 permeability.

It should be noted that since the velocity computations for thick layers are accurate for Da as low as 0.001, so are the computations of the mean velocity.

9. Velocity Profiles

Velocity profiles across the three-layered channel have been obtained for the various flow parameters and are graphed in Figures 3-9.

Figure 3 illustrates the velocity profile for a thick transition layer,

$\eta =\frac{1}{3},\xi =\frac{2}{3}$ , when Da = 1 and n = 1. This velocity profile is parabolic and is

close to what one might expect in the case of Poiseuille flow and the solution of the Navier-Stokes equations due to the high value of permeability. This parabolic shape is lost when Da is reduced, as shown in Figure 4, as the flow in upper part of the channel moves slower than in the middle part, which is in turn slower than the Navier-Stokes flow in the lower part of the configuration.

For the case of n = 1, and thick transition layer, Figure 5 illustrates the effect of changing Da on the velocity profile. The case of thin layer is shown in Figure 6. Both figures demonstrate a velocity profile distortion with decreasing Da as the permeability becomes lower in the upper region, resulting in slower flow. The same is illustrated in Figure 7, Figure 8 for the case of n = 2.

The effect of varying n is illustrated in Figure 9, which demonstrates the relative increase of velocity in the transition layer with increasing n. This may be

Figure 3. Velocity Profile u(y), for Da = 1, $\eta =\frac{1}{3},\xi =\frac{2}{3}$ , and n = 1.

Figure 4. Velocity Profile u(y), for Da = 0.01, $\eta =\frac{1}{3},\xi =\frac{2}{3}$ , and n = 1.

Figure 5. Velocity profiles in the three layers for n = 1, $\eta =\frac{1}{3},\xi =\frac{2}{3}$ and different Da.

Figure 6. Velocity profiles in the three layers for n = 1, $\eta =0.49,\xi =0.51$ and different values of Darcy number.

Figure 7. Velocity profiles in the three layers for n = 2, $\eta =\frac{1}{3},\xi =\frac{2}{3}$ and different Da.

Figure 8. Velocity profiles in the three layers for n = 2, $\eta =0.49,\xi =0.51$ and different Da.

Figure 9. Velocity profiles in the layers for Da = 0.1, $\eta =\frac{1}{3},\xi =\frac{2}{3}$ and different n.

attributed to higher momentum transfer from the adjacent layers that tend to increase the flow velocity there.

10. Conclusion

In this work, we considered the problem of flow through a porous medium over a free-space channel in the presence of a transition layer. This problem was treated by Nield and Kuznetsov [1] 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 considered the permeability to vary as the reciprocal of an nth degree polynomial and to solve the resulting generalized Airy’s differential equation. For the case of n = 1, our results agree with those obtained by Nield and Kuznetsov [1] .

Cite this paper

Zaytoon, M. , Alderson, T. and Hamdan, M. (2018) A Study of Flow through a Channel Bounded by a Brinkman Transition Porous Layer.*Journal of Applied Mathematics and Physics*, **6**, 264-282. doi: 10.4236/jamp.2018.61025.

Zaytoon, M. , Alderson, T. and Hamdan, M. (2018) A Study of Flow through a Channel Bounded by a Brinkman Transition Porous Layer.

References

[1] Nield, D.A. and Kuznetsov, A.V. (2009) The Effect of a Transition Layer between a Fluid and a Porous Medium: Shear Flow in a Channel. Transport in Porous Media, 78, 477-487.

https://doi.org/10.1007/s11242-009-9342-0

[2] Nield, D.A. (1983) The Boundary Correction for the Rayleigh-Darcy Problem: Limitations of the Brinkman Equation. Journal of Fluid Mechanics, 128, 37-46.

https://doi.org/10.1017/S0022112083000361

[3] Goyeau, B., Lhuillier, D. and Velarde, M.G. (2003) Momentum Transport at a Fluid-Porous Interface. International Journal of Heat and Mass Transfer, 46, 4071-4081.

https://doi.org/10.1016/S0017-9310(03)00241-2

[4] Goharzadeh, A., Khalili, A. and Jorgensen, B.B. (2005) Transition Layer Thickness at a Fluid-Porous Interface. Physics of Fluids, 17, 057102.

https://doi.org/10.1063/1.1894796

[5] Parvazinia, M., Nassehi, V., Wakeman, R.J. and Ghoreishy, M.H.R. (2006) Finite Element Modelling of Flow through a Porous Medium between Two Parallel Plates Using the Brinkman Equation. Transport in Porous Media, 63, 71-90.

https://doi.org/10.1007/s11242-005-2721-2

[6] Chandesris, M. and Jamet, D. (2007) Boundary Conditions at a Fluid-Porous Interface: An a priori Estimation of the Stress Jump Boundary Conditions. International Journal of Heat and Mass Transfer, 50, 3422-3436.

https://doi.org/10.1016/j.ijheatmasstransfer.2007.01.053

[7] Hill, A.A. and Straughan, B. (2008) Poiseuille Flow in a Fluid Overlying a Porous Medium. Journal of Fluid Mechanics, 603, 137-149.

https://doi.org/10.1017/S0022112008000852

[8] Duman, T. and Shavit, U. (2010) A Solution of the Laminar Flow for a Gradual Transition between Porous and Fluid Domains. Water Resources Research, 46, 09517.

https://doi.org/10.1029/2009WR008393

[9] Tao, J., Yao, J. and Huang, Z. (2013) Analysis of the Laminar Flow in a Transition Layer with Variable Permeability between a Free-Fluid and a Porous Medium. Acta Mechanica, 224, 1943-1955.

https://doi.org/10.1007/s00707-013-0852-z

[10] Hamdan, M.H. and Kamel, M.T. (2011) On the Ni(x) Integral Function and Its Application to the Airy’s Non-Homogeneous Equation. Applied Mathematics and Computation, 217, 7349-7360.

https://doi.org/10.1016/j.amc.2011.02.025

[11] Nield, D.A. and Bejan, A. (2006) Convection in Porous Media. 3rd Edition, Springer.

[12] Nield, D.A. (2009) The Beavers Joseph Condition and Related Matters: A Historical and Critical Note. Transport in Porous Media, 78, 537-540.

https://doi.org/10.1007/s11242-009-9344-y

[13] Swanson, C.A. and Headley, V.B. (1967) An Extension of Airy’s Equation. SIAM Journal on Applied Mathematics, 15, 1400-1412.

https://doi.org/10.1137/0115123

[1] Nield, D.A. and Kuznetsov, A.V. (2009) The Effect of a Transition Layer between a Fluid and a Porous Medium: Shear Flow in a Channel. Transport in Porous Media, 78, 477-487.

https://doi.org/10.1007/s11242-009-9342-0

[2] Nield, D.A. (1983) The Boundary Correction for the Rayleigh-Darcy Problem: Limitations of the Brinkman Equation. Journal of Fluid Mechanics, 128, 37-46.

https://doi.org/10.1017/S0022112083000361

[3] Goyeau, B., Lhuillier, D. and Velarde, M.G. (2003) Momentum Transport at a Fluid-Porous Interface. International Journal of Heat and Mass Transfer, 46, 4071-4081.

https://doi.org/10.1016/S0017-9310(03)00241-2

[4] Goharzadeh, A., Khalili, A. and Jorgensen, B.B. (2005) Transition Layer Thickness at a Fluid-Porous Interface. Physics of Fluids, 17, 057102.

https://doi.org/10.1063/1.1894796

[5] Parvazinia, M., Nassehi, V., Wakeman, R.J. and Ghoreishy, M.H.R. (2006) Finite Element Modelling of Flow through a Porous Medium between Two Parallel Plates Using the Brinkman Equation. Transport in Porous Media, 63, 71-90.

https://doi.org/10.1007/s11242-005-2721-2

[6] Chandesris, M. and Jamet, D. (2007) Boundary Conditions at a Fluid-Porous Interface: An a priori Estimation of the Stress Jump Boundary Conditions. International Journal of Heat and Mass Transfer, 50, 3422-3436.

https://doi.org/10.1016/j.ijheatmasstransfer.2007.01.053

[7] Hill, A.A. and Straughan, B. (2008) Poiseuille Flow in a Fluid Overlying a Porous Medium. Journal of Fluid Mechanics, 603, 137-149.

https://doi.org/10.1017/S0022112008000852

[8] Duman, T. and Shavit, U. (2010) A Solution of the Laminar Flow for a Gradual Transition between Porous and Fluid Domains. Water Resources Research, 46, 09517.

https://doi.org/10.1029/2009WR008393

[9] Tao, J., Yao, J. and Huang, Z. (2013) Analysis of the Laminar Flow in a Transition Layer with Variable Permeability between a Free-Fluid and a Porous Medium. Acta Mechanica, 224, 1943-1955.

https://doi.org/10.1007/s00707-013-0852-z

[10] Hamdan, M.H. and Kamel, M.T. (2011) On the Ni(x) Integral Function and Its Application to the Airy’s Non-Homogeneous Equation. Applied Mathematics and Computation, 217, 7349-7360.

https://doi.org/10.1016/j.amc.2011.02.025

[11] Nield, D.A. and Bejan, A. (2006) Convection in Porous Media. 3rd Edition, Springer.

[12] Nield, D.A. (2009) The Beavers Joseph Condition and Related Matters: A Historical and Critical Note. Transport in Porous Media, 78, 537-540.

https://doi.org/10.1007/s11242-009-9344-y

[13] Swanson, C.A. and Headley, V.B. (1967) An Extension of Airy’s Equation. SIAM Journal on Applied Mathematics, 15, 1400-1412.

https://doi.org/10.1137/0115123