Flow Governed by Generalized Brinkman’s Equation through an Inclined Porous Channel

1. Introduction

The main objective of this work is to study the effects of viscosity variations due to pressure on the flow through a porous structure of an incompressible fluid, governed by the generalized Brinkman equation. In particular, viscosity, $\mu$, is taken as an exponential function of pressure, $p$, of the form:

$$\mu \left(p\right)=f\left(p\right)={\mu }_0{e}^{\left({\left(\frac{p}{p_0}\right)}^n-1\right)}$$

(1)

where $\alpha \ge 0$, $p_0>0$, $n\ge -1$, and ${\mu }_0=f\left(p_0\right)>0$ is the constant viscosity when $\alpha =0$ and $p=p_0$.

The realization that changes in flow conditions (such as changes in temperature and pressure) result in corresponding changes in the viscosity of a flowing fluid dates back to the works of Stokes, [1], and Barus, [2]. With increasing temperature, viscosity of a fluid decreases, while an excessive increase in pressure could result in an increase in fluid viscosity. Barus, [2], suggested that viscosity is an exponentially increasing function of pressure of the form:

$$\mu ={\mu }_0{e}^{\alpha \left(p-p_0\right)}$$

(2)

For small values of $\alpha$ or small pressure differences, relationship (2) between viscosity and pressure can be expressed as, [3]:

$$\mu ={\mu }_0\left[1+\alpha \left(p-p_0\right)\right]$$

(3)

It has been customary to assume that in the flow of Newtonian liquids in the absence of temperature variations, shear viscosity is constant (cf. [4] and the references therein). This assumption ignores the effects of pressure on viscosity when the fluid pressure is of the order of 1000 atm., [3,5].

While the above relationships, (2) and (3), are suitable for fluids with small molecules, long chain molecules (the likes of polymers and some oil mixtures) require other forms of viscosity-pressure relations, [3,6,7]. It has also been realized that the effect of pressure on density of a fluid is small as compared to the effect of pressure on viscosity, [6,7,8]. Accordingly, it suffices in many applications to study incompressible fluid flow with pressure-dependent viscosity as governed by the equation of continuity and the Navier-Stokes equations, written here for steady flow, in the absence of body forces, in the following forms, [3,9,10], respectively:

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

(4)

$$\rho \overrightarrow{u}·\nabla \overrightarrow{u}=-\nabla p+\nabla ·\overrightarrow{T}+\rho \overrightarrow{g}$$

(5)

where

$$\overrightarrow{T}=-p\overrightarrow{I}+2\mu \left(p\right)\overrightarrow{K}$$

(6)

is the Cauchy stress tensor in which

$$\overrightarrow{K}=\frac{1}{2}(\nabla \overrightarrow{u}+\left(\nabla \overrightarrow{u}{)}^T\right)$$

(7)

is the symmetric part of the velocity gradient, $\overrightarrow{u}$ is the velocity field, $\rho$ is the constant fluid density, and $\rho \overrightarrow{g}$ is the gravitational force per unit volume.

Governing equations (4) and (5), with $\overrightarrow{T}$ and $\overrightarrow{K}$ defined by (6) and (7), represent an underdetermined system of four equations in the five unknowns $\overrightarrow{u}=\left(u,v,w\right)$, $p$ and $\mu$. In the absence of additional conservation principles to provide an additional equation, viscosity can be expressed as a function of pressure (such as in (2) and (3)) to provide an additional condition to render the governing system of equations determinate.

In addition to the flow conditions’ influence on viscosity, there are flow domain conditions that could result in changes in viscosity of a flowing fluid. For instance, in fluid flow through domains with constrictions and in porous structures, it has been suggested by a large number of researchers that the effective viscosity of a fluid in a porous structure is different than the viscosity of the base fluid (cf. [11] and the references therein). It can be argued that Brinkman’s equation accounts for viscosity change due to domain constrictions by introducing an effective viscosity term. However, it does not explicitly take into account spatial variations of viscosity, or viscosity variations due to changes in pressure, which is a function of position, (cf. [12] and the references therein).

When the flow is taken through a porous medium, Navier-Stokes equations are valid microscopically in the pore space. However, they are hard to track due to the complexity of the pore space. Furthermore, changes in flow quantities at the microscopic scale are at times insignificant at the macroscopic scale. The equations governing the flow of a fluid with pressure-dependent viscosity through a porous structure are the continuity equation, of the form of (4), and the momentum equation proposed by Rajagopal et.al, [10], referred to as the generalized Brinkman’s equation. For the purposes of this work, the generalized Brinkman’s equation is conveniently written in the following form:

$$\rho \left(\overrightarrow{u}·\nabla \right)\overrightarrow{u}=-\nabla p+\nabla ·\mu \left(\nabla \overrightarrow{u}+{\left(\nabla u\right)}^T\right)-\mathsf{\Gamma }\left(p\right)\overrightarrow{u}+\rho \overrightarrow{g}$$

(8)

where $\overrightarrow{u}$ in (4) and (8) now stands for the Brinkman velocity vector, $\mu =\mu \left(p\right)$, and $\mathsf{\Gamma }\left(p\right)$ is the pressure-dependent drag coefficient. Both $\mu \left(p\right)$ and $\mathsf{\Gamma }\left(p\right)$ control variations in viscosity due to pressure, and variations in pressure due to variations in porous parameters. In particular, Subramanian and Rajagopal, [13], identified $\mu \left(p\right)$ with the measure of resistance between fluid layers while $\mathsf{\Gamma }\left(p\right)$ is a measure of friction between solid and fluid at the pore.

Various forms and combinations of $\mu \left(p\right)$ and $\mathsf{\Gamma }\left(p\right)$ have been suggested and used in the literature, popular among which are the following forms, [14]:

$$\mu \left(p\right)=A{e}^{ap}\mathsf{\Gamma }\left(p\right)=B{e}^{bp}$$

(9)

$$\mu \left(p\right)=A{e}^{ap}\mathsf{\Gamma }\left(p\right)=B{\left(\frac{p}{p_0}\right)}^m$$

(10)

$$\mu \left(p\right)=A{\left(\frac{p}{p_0}\right)}^n\mathsf{\Gamma }\left(p\right)=B{e}^{bp}$$

(11)

$$\mu \left(p\right)=A{\left(\frac{p}{p_0}\right)}^n\mathsf{\Gamma }\left(p\right)=B{\left(\frac{p}{p_0}\right)}^m$$

(12)

wherein $A$, $B$, $a,b,m,n$ are constants.

Considering that flow through a porous structure is influenced in part by permeability to the flow, Alharbi et.al., [11], suggested that $\mathsf{\Gamma }\left(p\right)=\frac{\mu \left(p\right)}{k}$, where k is the permeability, and provided the following form of momentum equation (8):

$$\rho \left(\overrightarrow{u}·\nabla \right)\overrightarrow{u}=-\nabla p+\nabla ·\mu \left(\nabla \overrightarrow{u}+{\left(\nabla \overrightarrow{u}\right)}^T\right)-\frac{\mu \left(p\right)}{k}\overrightarrow{u}+\rho \overrightarrow{g}$$

(13)

It is conceivable that the use of (8) or (13) is predicated by the type of practical application being considered, and by the porous microstructure. In general, flow through porous media has well established applications that include the study of groundwater movement, oil and gas recovery, flow through membranes and kidney dialysis, air conditioning and heat transfer, and the drying of solids, (cf. [15] and the references therein). Flow with high pressures arise in a number of industrial applications that involve chemical and process technologies, such as medical tablet production, crude and fuel oil pumping, food processing, fluid-film lubrication theory, and in microfluidics, [3 ,4,5,16,17,18,19,20]. The form of governing equations used must accommodate and describe the different types of fluid flow, flow conditions, and the anticipated porous structure, [12,21,22].

In the current work, the objective is to study the effects of porous medium permeability, and the form of dependence of viscosity on pressure, on characteristics of the flow. In particular, relationship (1) between viscosity and pressure is employed in this work in order to study the effects of $n$. Cases where $n$ is a positive integer are investigated, together with the case when $n=-1$, which results in a decrease in the exponent with increasing pressure distribution. For comparison, the case of $n=0$ is considered. This case corresponds to flow with constant fluid viscosity.

To accomplish the above objective, model equation (13) will be used under the assumption of a constant permeability at this stage, while the more realistic variable permeability considerations will be considered in future work and after quantification of constant permeability influence on the flow. In order to quantify the effects of relationship (1) on the flow characteristics, unidirectional flow through an inclined porous channel is considered under variations in six parameters. The six parameters to be investigated are the angle of inclination of the channel to the horizontal; pressure-control parameter, $\alpha$; medium permeability, $k$; constant reference pressure, $p_0$ (taken as atmospheric pressure at upper channel wall); constant reference viscosity, ${\mu }_0$, taken as the value of viscosity at pressure $p_0$; and variations in $n$ (namely, the values -1, 0, 1, 2, and 3).

All solutions and simulations used in this work have been carried out using MATLAB software package, version R2022a. However, for the sake of comparison and validation of the numerical results, the analytic solution to the flow problem at hand is obtained for the case of $n=1$. Comparisons of flow characteristics obtained using the analytic solution and those obtained using MATLAB are made.

2. Problem Formulation

Consider the flow of a fluid with pressure-dependent viscosity through a porous sediment of depth h inclined at angle $\vartheta$ to the horizontal. The flow configuration is illustrated in Figure 1 which shows the orientation of the coordinate system used. It is assumed that the sediment is bounded by impermeable, solid walls on which the no-slip condition is applied.

Flow in the above domain is governed by the equation of continuity (4) and momentum equations (13), which reduce to the following set of equations, wherein $u=u\left(y\right)$ is the tangential velocity component and $g$ is acting vertically downward:

$$-\frac{dp}{dx}+\mu \frac{d^{2}u}{{dy}^2}+\frac{d\mu }{dy}\frac{du}{dy}+\rho gsin\vartheta -\frac{\mu }{k}u=0$$

(14)

$$-\frac{dp}{dy}-\rho gcos\vartheta =0$$

(15)

with boundary conditions given by zero-slip on the solid boundaries $y=0$ and $y=h$, and a prescribed pressure at $y=h$ (such as atmospheric pressure, say$p_0$). Boundary conditions are thus given as:

$$u\left(0\right)=0,u\left(h\right)=0,p\left(h\right)=p_0$$

(16)

Introducing the dimensionless quantities defined as:

$$y^{*}=\frac{y}{h};u^{*}=\frac{u}{\sqrt{gh}};p^{*}=\frac{p}{\rho gh};{\mu }^{*}=\frac{\mu }{\rho h\sqrt{gh}};k^{*}=\frac{k}{h^2}$$

(17)

then boundary conditions, (16), and governing equations (14) and (15) take the following dimensionless forms, respectively, after dropping the asterisks “*”:

$$u\left(0\right)=0,u\left(1\right)=0,p\left(1\right)=p_0$$

(18)

$$\mu \frac{d^{2}u}{{dy}^2}+\frac{d\mu }{dy}\frac{du}{dy}+sin\vartheta -\frac{\mu }{k}u=0$$

(19)

$$\frac{dp}{dy}=-cos\vartheta$$

(20)

We note that in (17), the velocity scale has been set to $\sqrt{gh}$, which renders the resulting solution readily applicable as the velocity scale depends on physical parameters of the problem. Furthermore, since $\sqrt{gh}$ is the velocity scale, $u^{*}$ is the local Froude number and ${\mu }^{*}=\frac{1}{Re}$, where $Re$ is the Reynolds number.

The general solution to (20) takes the form

$$p=c-ycos\vartheta$$

(21)

where c is an arbitrary constant.

Using pressure condition $p\left(1\right)=p_0$ we find that $c=p_0+cos\vartheta$ and (21) takes the form

$$p=p_0+\left(1-y\right)cos\vartheta$$

(22)

In order to solve (19) for $u\left(y\right)$, we assume that the viscosity varies with pressure. We thus let

$$\mu \left(p\right)=f\left(p\right)$$

(23)

where $f\left(p\right)$ is a function of $p$ (to be specified).

From (23), we obtain

$$\frac{d\mu }{dy}=\frac{d\mu }{dp}\frac{dp}{dy}=f^{\prime }\left(p\right)\frac{dp}{dy}$$

(24)

From (22), we obtain

$$\frac{dp}{dy}=-cos\vartheta$$

(25)

Using (25) in (24), we obtain

$$\frac{d\mu }{dy}={-cos\vartheta f}^{\prime }\left(p\right)$$

(26)

Using (23) and (26) in (19), we obtain

$$f\left(p\right)\frac{d^{2}u}{{dy}^2}{-cos\vartheta f}^{\prime }\left(p\right)\frac{du}{dy}+sin\vartheta -\frac{f\left(p\right)}{k}u=0$$

(27)

or

$$\frac{d^{2}u}{{dy}^2}-cos\vartheta \frac{f^{\prime }\left(p\right)}{f\left(p\right)}\frac{du}{dy}-\frac{u}{k}=-\frac{sin\vartheta }{f\left(p\right)}$$

(28)

Equation (28) is a general differential equation that governs the flow of a fluid with pressure-dependent viscosity through a porous domain down an inclined plane. Given $f\left(p\right)$ and $k$, specific forms of (28) are obtained. A form of interest to this work from a fundamental analysis point of view is the following:

$$f\left(p\right)={\mu }_0{e}^{\alpha \left({\left(\frac{p}{p_0}\right)}^n-1\right)}={A\mu }_0{e}^{{\left(\frac{p}{p_0}\right)}^n}$$

(29)

where $A=e^{-\alpha }$, $\alpha >0,{\mu }_0=f\left(p_0\right)>0\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y}\mathrm{a}\mathrm{t}\mathrm{y}=1,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}p=p_0>0$, and $n\ge -1$ is an integer.

Using (29) in (28), the following governing differential equation is obtained:

$$\frac{d^{2}u}{{dy}^2}-cos\vartheta \frac{\alpha n(p{)}^{n-1}}{{(p}_0{)}^n}\frac{du}{dy}-\frac{u}{k}=-\frac{sin\vartheta }{A{\mu }_0{e}^{\alpha {\left(\frac{p}{p_0}\right)}^n}}$$

(30)

Letting

$$P=\frac{p}{p_0}$$

(31)

then (17) can be written as:

$$P=\frac{p}{p_0}=B-Cy$$

(32)

where $B=1+C$ and $C=\frac{cos\vartheta }{p_0}$.

Using (31) and (32), equation (30) can be written as:

$$\frac{d^{2}u}{{dy}^2}-\alpha nC(B-Cy{)}^{n-1}\frac{du}{dy}-\frac{u}{k}=-\frac{sin\vartheta }{A{\mu }_0{e}^{\alpha (B-Cy{)}^n}}$$

(33)

Equation (33) is to be solved subject to no-slip conditions in (16). Solution to (33) is obtained numerically using MATLAB. However, when $n=1$, (33) reduces to

$$\frac{d^{2}u}{{dy}^2}-\alpha C\frac{du}{dy}-\frac{u}{k}=-\frac{sin\vartheta }{A{\mu }_0{e}^{\alpha B}}{e}^{\alpha Cy}$$

(34)

whose general solution takes the form:

$$u=b_1{e}^{m_{1}y}+b_2{e}^{m_{2}y}+\frac{ksin\vartheta e^{\alpha Cy}}{A{\mu }_0{e}^{\alpha B}}$$

(35)

Using (35) and (16), the following solution to the boundary value problem is obtained:

$$u=\frac{ksin\vartheta }{A{\mu }_0{e}^{\alpha B}}\left\{\left[\frac{e^{\alpha C}{-e}^{m_2}}{e^{m_2}-e^{m_1}}\right]e^{m_{1}y}+\left[\frac{{e^{m_1}-e}^{\alpha C}}{e^{m_2}-e^{m_1}}\right]e^{m_{2}y}+e^{\alpha Cy}\right\}$$

(36)

where

$$m_1=\frac{\alpha C-\sqrt{{\alpha }^2{C}^2+\frac{4}{k}}}{2}$$

(37)

$$m_2=\frac{\alpha C+\sqrt{{\alpha }^2{C}^2+\frac{4}{k}}}{2}$$

(38)

Using (36), vorticity, $\omega$, is obtained as:

$$\omega =-\frac{du}{dy}=-\frac{m_{1}ksin\vartheta }{A{\mu }_0{e}^{\alpha B}}\left\{\left[\frac{e^{\alpha C}{-e}^{m_2}}{e^{m_2}-e^{m_1}}\right]e^{m_{1}y}-m_2\left[\frac{{e^{m_1}-e}^{\alpha C}}{e^{m_2}-e^{m_1}}\right]e^{m_{2}y}-\alpha C{e}^{\alpha Cy}\right\}$$

(39)

and the shear stress, $\tau$, takes the following form:

$$\tau =\mu \frac{du}{dy}=-\mu \omega ={-A\mu }_0{e}^{\alpha (B-Cy{)}^n}\omega$$

(40)

Total shear stress is given by

$$\frac{d}{dy}\left(\mu \frac{du}{dy}\right)=\frac{d}{dy}(-\mu \omega )={\frac{d}{dy}\{-A\mu }_0{e}^{\alpha (B-Cy{)}^n}\omega \}$$

(41)

3. Results and Discussion

The solution given by equation (36) was tabulated for the Reference Case given in Table 1 and compared to the numerical solution calculated with MATLAB. The magnitude of the difference in the values of $u$ between the two solutions across the channel was less than $4\times 10^{-17}$.

Equation (33) was solved numerically using the bvp5c algorithm in MATLAB R2022a which is described by Kierzenka and Shampine [23] with the relative and absolute tolerances set to $1\times {10}^{-12}$. Results have been obtained for the ranges of flow and domain parameters given in Table 1:

Table 1. Ranges of Flow and Domain Parameters.

Case	1	2	3	4	5	6
Parameter			Reference Case
$\mathit{\vartheta }$	${15}^{°}$	${30}^{°}$	${45}^{°}$	${60}^{°}$	${75}^{°}$	${90}^{°}$
${\mathit{\mu }}_0$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
${\mathit{p}}_0$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
$\mathit{\alpha }$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
$\mathit{n}$	-1	0	1	2	3
$\mathit{k}$	${10}^{-4}$	${10}^{-3}$	${10}^{-2}$	${10}^{-1}$	1	-

Table 1. Ranges of Flow and Domain Parameters.

Case	1	2	3	4	5	6
Parameter			Reference Case
$\mathit{\vartheta }$	${15}^{°}$	${30}^{°}$	${45}^{°}$	${60}^{°}$	${75}^{°}$	${90}^{°}$
${\mathit{\mu }}_0$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
${\mathit{p}}_0$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
$\mathit{\alpha }$	${10}^{-1}$	${10}^{-0.5}$	1	${10}^{0.5}$	${10}^1$	-
$\mathit{n}$	-1	0	1	2	3
$\mathit{k}$	${10}^{-4}$	${10}^{-3}$	${10}^{-2}$	${10}^{-1}$	1	-

The column with Reference Case heading in Table 1 is used to isolate phenomena and contains the fixed values of parameters when a particular parameter is under investigation. The graphical results for the Reference Case appear in all of the parametric studies to facilitate comparison. For instance, to study the effects of changes in $\alpha$, its values are varied as listed in the horizontal row labelled $\alpha$ in Table 1, namely $\alpha$ takes the range of values 10^-1, 10^-0.5, 1, 10^0.5, and 10¹, while other parameters’ fixed values are listed in the reference case column, namely $\vartheta ={45}^{°}$, ${\mu }_0=1$, $p_0=1$, $n=1$, and $k={10}^{-2}$.

In the discussion that follows, effects of the parameters of Table 1 on viscosity and pressure distributions, velocity profiles, shear stress and viscous drag and resistance to the flow are discussed.

Six parametric studies will be presented in the following order: $\vartheta ,{\mu }_0,p_0,\alpha ,n$, and $k$. The results for each parametric study include (1) velocity profiles, $u\left(y\right)$; (2) viscosity profiles, $\mu \left(y\right)$; (3) Darcy friction profiles, $-\mu \left(y\right)\left(\frac{u\left(y\right)}{k}\right)$; (4) the net force on a differential fluid element due to variation of the shear strain rate across the fluid element, $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$; (5) the net force on a differential fluid element due to variation of the viscosity across the fluid element, $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$; and (6) the dimensionless volumetric flowrate in the channel, $Q$.

Items (3)-(5) show the spatial variation each of the forces present in the dimensionless force balance given by equation (19) and must sum to $-\sin \left(\vartheta \right)$, the component of the fluid’s weight acting along the incline. Presentation of the results in this manner aids in the development of an understanding of the physical phenomena involved in the resulting velocity profiles. Note that in some instances semi-logarithmic axes are used to better illustrate the results when the quantity in question varies over orders of magnitude.

A visual display of the ranges of Table 1 is shown in Figure 2, below, to demonstrate increase and decrease of the chosen values for the parameters relative to each other and how they influence the dimensionless volumetric flowrate in the channel, $Q$. Figure 2 shows that an increases in ${\mu }_0$ results in a linear decrease in $Q$. It also illustrates the expected largest rate of increase in $Q$ with increasing permeability, $k$.

3.1. Pressure and Viscosity Distributions

Prior to discussion of the resulting velocity and force profiles, the effects of the parameters on the pressure distribution, and hence the viscosity distributions, will be presented. Plots of the pressure across the channel are not included as they vary linearly according to equation (22).

In this work viscosity is taken to be a function of pressure as defined by equation (29), wherein parameters $n$ and $\alpha$ are the viscosity distribution control parameters. Therefore, as the pressure varies due to hydrostatic forces, the fluid’s viscosity also varies across the channel. In practice, the values of $p_0$ and ${\mu }_0$, $n$ and $\alpha$ can adjusted to reproduce the appropriate variation of the viscosity for the problem at hand. In the present study, their chosen ranges guarantee realistic viscosity values yet are varied enough to illustrate their effects on the flow. Other than the permeability, $k$, variations in all other parameters listed in Table 1 have a direct influence on viscosity across the channel, as discussed in what follows.

3.2. Effects of $\vartheta$

Viscosity is assumed to be a function of pressure as given by (29) while the pressure in the channel is given by (22), namely $p\left(y\right)=p_0+\left(1-y\right)cos\vartheta$. From this equation it is observed that the pressure has its smallest value, $p_0$, at the upper wall and increases linearly with decreasing $y$ until it reaches its maximum $p_0+cos\vartheta$ at the lower wall. Changes in pressure across the channel affect the viscosity distribution across the channel, and result in a maximum viscosity at the lower wall and a minimum at the upper wall (in all but one case, when $n=-1$). Effects of $\vartheta$ on $u,$ $\mu$, $Q$, $-\frac{\mu u}{k}$, $\mu \frac{d^{2}u}{d{y}^2}$, and $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$ are illustrated in Figure 3.

The inclinations of the flow domain used were $\vartheta =15°,30°,45°,60°,75°,\mathrm{a}\mathrm{n}\mathrm{d}90°$ while all other parameters take on the values identified in the Reference Case given in Table 1. The variation of the angle has two distinct effects on the flow: the first is the variation of the pressure across the channel, and thereby the variation of the viscosity; the second effect is the variation of the component of gravitational acceleration along the incline, and therefor the driving force for the flow. When $\vartheta =0°$ the variation of the pressure across the channel is maximized but does not produce flow as there is no component of gravity acting along the incline. When the angle of inclination is ${90}^{°}$, equations (22) and (29) render $p\left(y\right)=p_0$ and $\mu \left(y\right)={\mu }_0$ meaning the viscosity is constant across the channel.

The plots of the velocity profiles, which can be found in Figure 3 show that for the case of $\vartheta =90°$, the velocity profile is symmetric about the center of the channel. This is a result of the unform pressure, and therefore uniform viscosity, across the channel. As $\vartheta$ decreases, the component of gravitational acceleration along the direction of flow decreases thereby reducing the overall driving force behind the flow. Furthermore, the higher viscosity as at the bottom of the channel, illustrated in Figure 3 causes the velocity to be lower near the bottom of the channel resulting in the asymmetry in the velocity profiles which becomes more pronounced at smaller values of $\vartheta$. Despite the lower velocity, Figure 3 shows that the Darcy friction is significant and nearly symmetric across the channel with the low-speed regions corresponding to regions of high viscosity. As mentioned previously, when $\vartheta =90°$ there is no variation in pressure, and therefore viscosity, across the channel so $\frac{d\mu }{dy}=0$ and a symmetric velocity profile and, therefor, symmetric distributions of Darcy friction and $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$. Examining the cases of $\vartheta \ne 90°$, small asymmetries can be seen in all of the forces and in the case of $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$, changes in direction of the force can be seen through the channel.

3.3. Effects of ${\mu }_0$

Equation (29) suggests that ${\mu }_0$ is a “multiplier” of $\mu \left(y\right)$, the dimensionless viscosity, across the channel. At $y=1$, or the upper channel wall, $\mu \left(1\right)={\mu }_0$, therefore, increasing ${\mu }_0$ results in increasing $\mu \left(1\right)$. At the lower channel wall, $y=0$, values of $\mu \left(0\right)$ are also dependent on ${\mu }_0$, as calculated using (29) and (32), for a given value of $\vartheta$.

The results for logarithmically evenly spaced values of ${\mu }_0$ can be found in Figure 4 where ${\mu }_0=10^{-1},{10}^{-0.5},1,{10}^{0.5},\mathrm{a}\mathrm{n}\mathrm{d}10^1$ while all other parameters take on the values identified in the Reference Case given in Table 1. Note that the velocity profiles are shown on a semi-log scale due to the significant reduction caused by the high viscosity flows. In terms of $\mu$ being the inverse of Reynolds number, an increase in $\mu$ implies a decrease in Reynolds number, hence a slower flow. As seen in Figure 4, the variation of ${\mu }_0$ produces a scaling of the velocity profile (and therefore $(\frac{du}{dy}$) and $\left(\frac{d^{2}u}{d{y}^2}\right)$) by a factor of $\frac{1}{{\mu }_0}$. This is a result of each of these quantities being multiplied by ${\mu }_0$ in equation (27) or, as seen more directly in equation (33), the division of the RHS by ${\mu }_0$. These effects combine to produce identical Darcy friction and viscous forces profiles for all values of ${\mu }_0$. As a result of the velocity profiles scaling by the factor of $\frac{1}{{\mu }_0}$, the dimensionless flowrate, Q, becomes a straight line when plotted on log-log axis. Since the profile of $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$ is not symmetric, the Darcy friction and $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$ also contain some slight asymmetry.

3.4. Effects of $p_0$ and $\alpha$

The value of $p_0$ is the value of pressure at the upper channel wall but also affects the pressure distribution across the channel as seen in equation (22). In equation (29), $p_0$ appears in the denominator of the exponent that defines viscosity.

Parameter $\alpha$ is a viscosity-control parameter which, when varied on its own, adjusts the rate at which the viscosity changes with pressure. As seen in equation (29), an increase in $\alpha$ results in an increase in the exponent that defines viscosity as a function of pressure. Since the value of the viscosity on the upper channel wall is set via ${\mu }_0$, increasing the value of $\alpha$ results in increasing the viscosity in the lower portion of the channel.

The effects of increasing $p_0$ on viscosity distribution are, therefor, opposite to the effects of increasing $\alpha$.

The results for the variation of $p_0$ and $\alpha$ can be seen in Figure 5 and Figure 6, respectively, while all other parameters take on the values identified in the Reference Case given in Table 1. Increasing $p_0$ translates into a higher value of pressure on the upper channel wall, due to pressure boundary condition there. At the lower channel wall, $y=0$, equation (22) gives ${p\left(0\right)=p}_0+cos\vartheta$ which, for a given $\vartheta$, increases with increasing $p_0$. Pressure distribution across the channel, therefore, decreases from ${p\left(0\right)=p}_0+cos\vartheta$ to ${p\left(1\right)=p}_0$. Clearly, the higher the value of $p_0$, the lower the ratio $P=\frac{p}{p_0}$ that appears in equation (29). This implies that the viscosity across the channel decreases and offers less resistnce to the flow. The end result is that flow velocity increases with increasing $p_0$.

The effect of low $p_0$ (high $\alpha$) is to cause the viscosity to vary more rapidly throughout the channel thereby decreasing the flow’s velocity in the lower, highly viscous, regions of the domain. Furthermore, this rapid variation of viscosity across the channel causes significant variation of the shear force across a differential fluid element as seen in the plot of $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$ in Figure 5, and a corresponding change in $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$; the latter developing a significant region of reversal of direction caused by the region of inflection in the velocity profile.

As $\alpha \to 0^{+}$, $\mu \left(y\right)\to {\mu }_0$ as can be seen from equation (29). When viscosity is constant across the channel, the velocity profile is near symmetric about the longitudinal centreline of the channel. As $\alpha$ increases, viscosity distribution across the channel decreases with increasing $y$, and the flow experiences more resistance in the lower parts of the channel. Velocity therefore decreases with decreasing $y$, accompanied with greater deviations from symmetry in the velocity profile. This behaviour is demonstrated in Figure 6 when $\alpha$ increases from 0.01 to 10.

3.5. Effects of $n$

Effects of variations in $n$ on the viscosity distribution $\mu \left(y\right)$ across the porous channel are illustrated in Figure 2 for the values of parameters identified in Table 1 as Reference Case. Figure 2 shows that at the upper channel wall, $y=1$, viscosity is independent of $n$. This can be seen from equation (29) where at the upper wall $p=p_0$, hence $\mu \left(1\right)=f\left(p_0\right)={\mu }_0=1$ for all values of $n$.

At the lower channel bounding wall, $y=0$, equation (22) gives the pressure $p\left(0\right)=p_0+cos{45}^{°}=1+\frac{1}{\sqrt{2}}.$ The quantity $\frac{p\left(0\right)}{p_0}=1+\frac{1}{\sqrt{2}}$results in increasing $\mu \left(y\right)$ with increasing $n$, as can be seen from equation (29). This translates into increasing $\mu \left(0\right)$ and an associated nonlinear increase in the viscosity profiles across the channel with increasing $n$.

The results from the variation of $n$ can be found in Figure 7 for $n=-1,0,1,2,\mathrm{a}\mathrm{n}\mathrm{d}3$ while all other parameters take on the values identified in the Reference Case given in Table 1. for the case of $n=0$, $\left(\frac{d\mu }{dy}\right)=0$ so the resulting flow is simply that of a constant viscosity fluid and corresponding symmetric velocity, Darcy friction, and $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$ profiles.

For the cases of $n>0$, $\mu$ increases with depth thereby producing a reduction in the velocity in the lower part of the channel. In the extreme case of $n=3$, a significant region of inflection is present in the velocity profile with a corresponding sign reversal of the viscous shear term $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$. Furthermore, the dramatic variation of the viscosity is also seen to produce significant variation of $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$ across the domain.

For the case of $n=-1$, $\left(\frac{d\mu }{dy}\right)>0$ meaning that, unlike all of the other flows considered in this work, the viscosity decreases with depth. The reduction in viscosity with depth results in less resistance to flow in the lower portion of the channel and, therefore, higher velocities there.

3.6. Variation of $k$

The results for the variation of the dimensionless permeability, $k$, can be found in Figure 8 while all other parameters take on the values identified in the Reference Case given in Table 1. Note that the velocity profiles are shown on a semi-log scale due to the significant reduction caused by the low permeability flows.

For the higher permeability flows the velocity profiles tend towards Navier-Stokes behaviour with a variable viscosity. As $k$ decreases, however, the resistance to flow through the porous medium increases resulting in a flow that is dominated by the Darcy friction. Taking equation (28) in the limit of small $k$ and then substituting for $f\left(p\right)$ from equation (29) and rearranging gives:

$$u\left(y\right)=\frac{k\sin \vartheta }{{A\mu }_0{e}^{{\left(\frac{p}{p_0}\right)}^n}}$$

(43)

which, for the Reference Case with $n=1$, will produce an inverse exponentially varying velocity profile across the channel. Near the edges of the channel, where the solution must satisfy the boundary conditions, the viscous shear terms become significant. As $k$ becomes progressively small, the region of flow dominated by Darcy-type behaviour increases. As seen in Figure 8 the Darcy friction profiles are relatively symmetric and, as $k$ becomes small, become nearly constant across the channel except near the boundaries. For small values of $k$ the viscous shear term, $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$, is nearly zero except near the walls and shows significant deviation from the constant value that would be expected in non-porous channel flow that would produce a parabolic velocity profile. The plot of the variation of $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$ can be found in Figure 8 where, since the variation of $\mu$ is the same in all flows, the differences in the profiles are due directly to the changes in the velocity profiles that result from the variation of $k$. Note that this force changes sign approximately at the midpoint of the channel where $\left(\frac{du}{dy}\right)$ changes sign to meet the no-slip boundary condition.

4. Conclusion

In this work, unidirectional flow of a fluid with pressure-dependent viscosity through a porous structure has been considered when the viscosity-pressure relationship is an exponential function of a pressure power function in order to investigate effects of the viscosity-pressure relation on the flow characteristics. This form of dependence is given by equation (1), above, when n takes on the values -1, 0, 1, 2, and 3. To accomplish this study, momentum equations governing the flow are given by (13), which have the realistic advantage of taking into account medium permeability in the drag term. It is realistic since one is dealing with flow through a porous structure. All solutions and simulations used in this work have been carried out using MATLAB software package, version R2022a. Effects of the six flow and medium parameters, $\vartheta ,{\mu }_0,p_0,\alpha ,n$, and $k$ on: (1) velocity profiles, $u\left(y\right)$; (2) viscosity profiles, $\mu \left(y\right)$; (3) Darcy friction profiles, $-\mu \left(y\right)\left(\frac{u\left(y\right)}{k}\right)$; (4) the net force on a differential fluid element due to variation of the shear strain rate across the fluid element, $\mu \left(\frac{d^{2}u}{d{y}^2}\right)$; (5) the net force on a differential fluid element due to variation of the viscosity across the fluid element, $\left(\frac{d\mu }{dy}\right)\left(\frac{du}{dy}\right)$; and (6) the dimensionless volumetric flowrate in the channel, $Q$, have been studied and analysed. Future work includes extension of the current constant permeability model to one of variable permeability, and contrasting the results with current ones.