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,
, is taken as an exponential function of pressure,
, of the form:
where
,
,
, and
is the constant viscosity when
and
.
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:
For small values of
or small pressure differences, relationship (2) between viscosity and pressure can be expressed as, [
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:
where
is the Cauchy stress tensor in which
is the symmetric part of the velocity gradient,
is the velocity field,
is the constant fluid density, and
is the gravitational force per unit volume.
Governing equations (4) and (5), with and defined by (6) and (7), represent an underdetermined system of four equations in the five unknowns , and . 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:
where
in (4) and (8) now stands for the Brinkman velocity vector,
, and
is the pressure-dependent drag coefficient. Both
and
control variations in viscosity due to pressure, and variations in pressure due to variations in porous parameters. In particular, Subramanian and Rajagopal, [
13], identified
with the measure of resistance between fluid layers while
is a measure of friction between solid and fluid at the pore.
Various forms and combinations of
and
have been suggested and used in the literature, popular among which are the following forms, [
14]:
wherein
,
,
are constants.
Considering that flow through a porous structure is influenced in part by permeability to the flow, Alharbi
et.al., [
11], suggested that
, where
k is the permeability, and provided the following form of momentum equation (8):
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 . Cases where is a positive integer are investigated, together with the case when , which results in a decrease in the exponent with increasing pressure distribution. For comparison, the case of 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, ; medium permeability, ; constant reference pressure, (taken as atmospheric pressure at upper channel wall); constant reference viscosity, , taken as the value of viscosity at pressure ; and variations in (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 . 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
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
is the tangential velocity component and
is acting vertically downward:
with boundary conditions given by zero-slip on the solid boundaries
and
, and a prescribed pressure at
(such as atmospheric pressure, say
). Boundary conditions are thus given as:
Introducing the dimensionless quantities defined as:
then boundary conditions, (16), and governing equations (14) and (15) take the following dimensionless forms, respectively, after dropping the asterisks “*”:
We note that in (17), the velocity scale has been set to , which renders the resulting solution readily applicable as the velocity scale depends on physical parameters of the problem. Furthermore, since is the velocity scale, is the local Froude number and , where is the Reynolds number.
The general solution to (20) takes the form
where
c is an arbitrary constant.
Using pressure condition
we find that
and (21) takes the form
In order to solve (19) for
, we assume that the viscosity varies with pressure. We thus let
where
is a function of
(to be specified).
Using (25) in (24), we obtain
Using (23) and (26) in (19), we obtain
or
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
and
, specific forms of (28) are obtained. A form of interest to this work from a fundamental analysis point of view is the following:
where
,
, and
is an integer.
Using (29) in (28), the following governing differential equation is obtained:
Letting
then (17) can be written as:
where
and
.
Using (31) and (32), equation (30) can be written as:
Equation (33) is to be solved subject to no-slip conditions in (16). Solution to (33) is obtained numerically using
MATLAB. However, when
, (33) reduces to
whose general solution takes the form:
Using (35) and (16), the following solution to the boundary value problem is obtained:
where
Using (36), vorticity,
, is obtained as:
and the shear stress,
, takes the following form:
Total shear stress is given by
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
between the two solutions across the channel was less than
.
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
. Results have been obtained for the ranges of flow and domain parameters given in
Table 1:
Table 1.
Ranges of Flow and Domain Parameters.
Table 1.
Ranges of Flow and Domain Parameters.
Case |
1 |
2 |
3 |
4 |
5 |
6 |
Parameter |
|
|
Reference Case
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1 |
|
|
- |
|
|
|
1 |
|
|
- |
|
|
|
1 |
|
|
- |
|
-1 |
0 |
1 |
2 |
3 |
|
|
|
|
|
|
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
, its values are varied as listed in the horizontal row labelled
in
Table 1, namely
takes the range of values 10
-1, 10
-0.5, 1, 10
0.5, and 10
1, while other parameters’ fixed values are listed in the reference case column, namely
,
,
,
, and
.
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: , and . The results for each parametric study include (1) velocity profiles, ; (2) viscosity profiles, ; (3) Darcy friction profiles, ; (4) the net force on a differential fluid element due to variation of the shear strain rate across the fluid element, ; (5) the net force on a differential fluid element due to variation of the viscosity across the fluid element, ; and (6) the dimensionless volumetric flowrate in the channel, .
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 , 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,
.
Figure 2 shows that an increases in
results in a linear decrease in
. It also illustrates the expected largest rate of increase in
with increasing permeability,
.
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
and
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
and
,
and
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,
, 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
Viscosity is assumed to be a function of pressure as given by (29) while the pressure in the channel is given by (22), namely
. From this equation it is observed that the pressure has its smallest value,
, at the upper wall and increases linearly with decreasing
until it reaches its maximum
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
). Effects of
on
,
,
,
, and
are illustrated in
Figure 3.
The inclinations of the flow domain used were
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
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
, equations (22) and (29) render
and
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
, 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
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
. 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
there is no variation in pressure, and therefore viscosity, across the channel so
and a symmetric velocity profile and, therefor, symmetric distributions of Darcy friction and
. Examining the cases of
, small asymmetries can be seen in all of the forces and in the case of
, changes in direction of the force can be seen through the channel.
3.3. Effects of
Equation (29) suggests that is a “multiplier” of , the dimensionless viscosity, across the channel. At , or the upper channel wall, , therefore, increasing results in increasing . At the lower channel wall, , values of are also dependent on , as calculated using (29) and (32), for a given value of .
The results for logarithmically evenly spaced values of
can be found in
Figure 4 where
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
being the inverse of Reynolds number, an increase in
implies a decrease in Reynolds number, hence a slower flow. As seen in
Figure 4, the variation of
produces a scaling of the velocity profile (and therefore
) and
) by a factor of
. This is a result of each of these quantities being multiplied by
in equation (27) or, as seen more directly in equation (33), the division of the RHS by
. These effects combine to produce identical Darcy friction and viscous forces profiles for all values of
. As a result of the velocity profiles scaling by the factor of
, the dimensionless flowrate,
Q, becomes a straight line when plotted on log-log axis. Since the profile of
is not symmetric, the Darcy friction and
also contain some slight asymmetry.
3.4. Effects of and
The value of 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), appears in the denominator of the exponent that defines viscosity.
Parameter 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 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 , increasing the value of results in increasing the viscosity in the lower portion of the channel.
The effects of increasing on viscosity distribution are, therefor, opposite to the effects of increasing .
The results for the variation of
and
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
translates into a higher value of pressure on the upper channel wall, due to pressure boundary condition there. At the lower channel wall,
, equation (22) gives
which, for a given
, increases with increasing
. Pressure distribution across the channel, therefore, decreases from
to
. Clearly, the higher the value of
, the lower the ratio
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
.
The effect of low
(high
) 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
in
Figure 5, and a corresponding change in
; the latter developing a significant region of reversal of direction caused by the region of inflection in the velocity profile.
As
,
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
increases, viscosity distribution across the channel decreases with increasing
, and the flow experiences more resistance in the lower parts of the channel. Velocity therefore decreases with decreasing
, accompanied with greater deviations from symmetry in the velocity profile. This behaviour is demonstrated in
Figure 6 when
increases from 0.01 to 10.
3.5. Effects of
Effects of variations in
on the viscosity distribution
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,
, viscosity is independent of
. This can be seen from equation (29) where at the upper wall
, hence
for all values of
.
At the lower channel bounding wall, , equation (22) gives the pressure The quantity results in increasing with increasing , as can be seen from equation (29). This translates into increasing and an associated nonlinear increase in the viscosity profiles across the channel with increasing .
The results from the variation of
can be found in
Figure 7 for
while all other parameters take on the values identified in the
Reference Case given in
Table 1. for the case of
,
so the resulting flow is simply that of a constant viscosity fluid and corresponding symmetric velocity, Darcy friction, and
profiles.
For the cases of , increases with depth thereby producing a reduction in the velocity in the lower part of the channel. In the extreme case of , a significant region of inflection is present in the velocity profile with a corresponding sign reversal of the viscous shear term . Furthermore, the dramatic variation of the viscosity is also seen to produce significant variation of across the domain.
For the case of , 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
The results for the variation of the dimensionless permeability,
, 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
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
and then substituting for
from equation (29) and rearranging gives:
which, for the
Reference Case with
, 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
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
becomes small, become nearly constant across the channel except near the boundaries. For small values of
the viscous shear term,
, 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
can be found in
Figure 8 where, since the variation of
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
. Note that this force changes sign approximately at the midpoint of the channel where
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, , and on: (1) velocity profiles, ; (2) viscosity profiles, ; (3) Darcy friction profiles, ; (4) the net force on a differential fluid element due to variation of the shear strain rate across the fluid element, ; (5) the net force on a differential fluid element due to variation of the viscosity across the fluid element, ; and (6) the dimensionless volumetric flowrate in the channel, , 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.
Figure 1.
Representative sketch.
Figure 1.
Representative sketch.
Figure 2.
Visual Display of Parameters’ Variations and their Effects on the dimensionless volumetric flowrate in the channel, .
Figure 2.
Visual Display of Parameters’ Variations and their Effects on the dimensionless volumetric flowrate in the channel, .
Figure 3.
Effects of on , , , , and .
Figure 3.
Effects of on , , , , and .
Figure 4.
Effects of on , , , , and .
Figure 4.
Effects of on , , , , and .
Figure 5.
Effects of on , , , , and .
Figure 5.
Effects of on , , , , and .
Figure 6.
Effects of on , , , , and .
Figure 6.
Effects of on , , , , and .
Figure 7.
Effects of on , , , , and .
Figure 7.
Effects of on , , , , and .
Figure 8.
Effects of on , , , , and .
Figure 8.
Effects of on , , , , and .