Finite Element Solutions for Non-Newtonian Pulsatile Flow in a Non-Darcian Porous Medium Conduit

Mathematics Department, Indian Institute of Technology Roorkee-247667, India rbharfma@iitr.ernet.in Engineering Department, Manchester Metropolitan University Oxford Rd., Manchester, M5 1GD, UK h.s.takhar@mmu.ac.uk Engineering Mechanics and Earthquake Dynamics Research 18 Milton Grove, Manchester, M16 OBP a.beg@email.com Fire Engineering Sciences, Leeds College of Building / Leeds Metropolitan University North Street, Leeds, LS2 7QT, UK obeg@lcb.ac.uk


Introduction
Pulsating flows abound in many areas of engineering fluid dynamics. Pressure surges in pipelines, cavitation in hydraulic systems [1], pumping of slurries and foodstuffs [2] etc are just several examples of the diverse applications of pulsatile hydromechanics. Other areas of technology in which pulsatile flows are important include refrigeration systems, combusting mechanisms, de-watering devices and also cardiovascular biomechanics. The reader is refereed to the excellent treatise by Zamir [3]. Pulsating flows are characterized by fluctuations in both mass flow rate and pressure about a non-zero mean value. Pulsatile flows can be generated by a variety of techniques. These include reciprocating pumps, steady flow pumps, hybrid systems with mechanical pulsation devices etc. An early study of pulsatile hydrodynamics was reported by Richardson and Tyler [4] who studied the alternating aerodynamic flow regimes in the vicinity of the entry region of pipes under periodic pressure differences. Sexl [5] later re-examined this problem and considered also the pulsatile flow generated by a reciprocating piston. Much later Uchida [6] considered oscillatory boundary layer flow regimes in pulsating pipe flows. These analyses provide classical benchmarks for more modern numerical studies, as documented by Schlichting [7]. Edwards and Wilkinson [8] later provided a succinct discussion of the applications of pulsatile flows. Other Newtonian studies of pulsating pipe flows include those by Ishii [9], Yakhot et al. [10] and most recently by Blyth and Pozrikidis [11] who studied numerically the influence of pulsatile flow on gas column stability. These studies were all concerned with purely fluid regimes. However it has been documented for some time that porous materials often occur in industrial operations and can be used to achieve hydrodynamic control of key processes. A lucid discussion of applications of porous media in conduit and channel flows is available in the excellent monograph by Scheidegger [12]. Wang [13] investigated the pulsatile Newtonian flow in a Darcian porous channel with suction/injection effects. In the present study we are concerned primarily with studying the flow of a non-Newtonian fluid under pulsatile pressure gradient. Several studies of pulsatile non-Newtonian flow through porous channels have appeared, mainly in the context of biomechanics. For example, Bhatnagar [14] analyzed theoretically the pulsating viscoelastic flow in a Darcian porous channel. The Darcian model [15] has been shown to be generally valid for low Reynolds number flows. At higher Reynolds numbers, inertial effects can become significant and the Darcian model is usually extended to incorporate a second order drag force. In the present study we therefore consider the pulsatile non-Newtonian flow in a Darcy-Forchheimer porous medium channel using a bi-viscosity rheological flow model described by Nakamura and Sawada [16]. The governing momentum equation is non-dimensionalized and solved subject to appropriate boundary conditions using the finite element technique. Several special cases are also discussed briefly. Such a study has to the author's knowledge not appeared thusfar in the technical literature. The present model may be applied to filtration processes in chemical engineering processes.

Hydromechanics of porous medium flow and rheological fluid model
Many techniques exist for simulating porous media. Adler [17] has described a wide range of geometric approaches which re-create the complex structure of porous materials.
In the context of engineering analysis we adopt instead a drag-force formulation for simulating the impedance imparted by the porous matrix structure to pulsating fluid in the channel. This approach has been applied in many porous media transport problems including geothermal systems, fibrous insulations, biomaterials and ceramics. Following Vafai and Tien [18] the appropriate formulation for flow via an isotropic, homogenous, non-stratified porous medium with high velocity effects i.e. the Darcy-Forchheimer model can be stated as follows: where V designates the velocity vector, p is the hydrodynamic pressure, µ is the dynamic viscosity, K is hydraulic conductivity (permeability) of the porous material, ρ is the density of the fluid and b is the Forchheimer geometrical (inertial drag) coefficient. A onedimensional version of this equation with velocity vector V = (U, 0, 0) will be employed in our model. The above equations relate to Newtonian porous flows. To accommodate the rheological characteristics of a slurry material, a stress tensor, τ ij , has to be introduced into the momentum equation. Following Nakamura and Sawada [16] we employ the biviscosity modified Casson rheological model. The constitutive equations for such a fluid may be presented as follows: where π = e ij e ij and e ij is the (i, j) th component of the deformation rate, π is the product of the component of deformation rate with itself, π c is a critical value of this product based on the Nakamura-Sawada model [16], µ B is plastic dynamic viscosity of the non-Newtonian fluid, p y is yield stress of slurry fluid. This model allows a relatively easy incorporation into the Navier-Stokes framework, as does the porous media formulation in (1).

Mathematical model
Consider the pulsatile, laminar, incompressible, flow of a Nakamura-Sawada fluid through a two-dimensional Darcy-Forchheimer porous medium channel possessing perforated walls. The physical regime is depicted in Fig. 1. The channel walls are located a distance H apart with reference to an x, y coordinate system, where x is in the longitudinal direction, parallel to the walls and y is normal to this. Pressure gradient is directed along the longitudinal axis. Transpiration can take place via the walls of the channel at a uniform velocity, V o . The convective terms in the momentum conservation equations i.e. u ∂u ∂x , u ∂v ∂x are ignored, following Uchida [6] which is physically allowable owing to the dominance of the viscous diffusion terms in pulsatile flow [3]. The y-direction velocity is neglected in the momentum equation and only simulated via a transpiration velocity where u is the longitudinal velocity component, t denotes time, ν is the kinematic viscosity of the slurry fluid, P is the pressure, ν B is plastic kinematic viscosity of fluid, β denotes the upper limit of apparent viscosity coefficient and all other parameters have been defined previously. The penultimate term in equation (3) For infinite values of β, the yield stress of the fluid vanishes and the flow becomes Newtonian. For the case of a channel with solid walls, V o → 0. We now introduce a series of dimensionless transformations to reduce the momentum equation with corresponding boundary conditions to non-dimensional form, aiding in a numerical solution to the flow problem. Defining: where X is the transformed coordinate parallel to channel walls, Y is the transformed coordinate transverse to channel walls, U is the transformed velocity component in x direction, V o is the transpiration velocity, P * is the transformed hydrodynamic pressure ( * dropped for convenience in analysis), T is dimensionless time. We therefore obtain the following second order partial differential equation for momentum in terms of dimensionless longitudinal velocity, U : where Re = HVo νB is the Reynolds number, λ = KVo νH is the Darcian permeability parameter and Nf = Hb is the Forchheimer inertial porous parameter. Effectively the transpiration has been absorbed into the Reynolds number and Darcian parameter and therefore no longer appears explicitly in the transformed model. The corresponding transformed boundary conditions now become: The pressure gradient in (6) is now decomposed into a steady component and an imposed (oscillatory) component as follows, as described by Uchida [6]:

Numerical solution
In the numerical (finite element) computations, the pulsating pressure gradient in the channel is re-defined as: Implementing this expression in (6) we arrive at: The initial condition in time is now prescribed as The entire flow domain is discretized into a set of 100 line elements of equal width, each element being two-noded. A variational formulation is employed in the finite element solution using a two-noded linear element with arbitrary test functions. Further details are provided in Bathe [19]. We are principally interested in studying velocity distributions in the channel in the spatial dimension (Y ) and over time (T ).
To study the interaction of pulsatile effects and rheological and non-Darcian effects in the flow regime we individually vary each of the following parameters: w * (dimensionless angular frequency), Re (Reynolds number), β (non-Newtonian parameter), λ (Darcian linear drag parameter), Po (pulsating amplitude), Ps (steady component of pressure gradient) and Nf (Forchheimer second order parameter). Velocity functions are plotted against Y and/or T (dimensionless time). Fig. 2 illustrates the visualization of velocity U with both space (Y coordinate) and time (T ). The oscillatory evolution of velocity is clearly visible. Computations correspond to a strongly non-Newtonian fluid. 3 peak velocities are observed within a time-span of 0 < T < 2. Peak velocities always correspond to the centerline of the channel i.e. at Y = 0.5. The undulating profile is a direct function of the parameters Po and Ps. Longitudinal velocity U starts from 0 at Y = 0 ascends to a peak at Y = 0.5 (centre of channel) and descends back to zero at Y = 1. The velocity also commences at 0 at T = 0 (initiation of the pulse) and the oscillatory flow nature is witnessed in time, rather than due to space, since the pulsatile gradient is a time-imposed effect. The computation corresponds to a strong pulsatile pressure gradient (w * = 8). Numerical solutions have been obtained for the general case by selecting a particular time T = 0.5 and computing the velocity function profiles at this instant. We utilize as default values, w * = 8 (angular frequency of oscillation), Re = 1 (Reynolds number), λ = 1 (Darcian linear drag parameter), Nf = 1 (Forchheimer second order parameter), β = 4 (non-Newtonian parameter), Ps = 10 (steady component of pressure gradient), Po = 7 (pulsating amplitude). All general cases correspond to these values. For special cases, corresponding values of physical parameters are provided in the following discussion, rather than on the graphs.  haviour. Increasing β decreases the effective viscosity of the fluid which tends to the Newtonian flow case. As β values tend to infinity, the term 1/β in equation (10) vanishes. Velocity U increases markedly with a rise in β; peak value of U at the centre of the channel rises from 15 to approximately 18 as β increases from 1 to 10 i.e. velocity is boosted by 20 %. The fluid flow tends further from shear-thickening behaviour as β is increased. The flow field is therefore clearly accelerated. In order to overcome the lower velocities for stronger rheological fluid, pulsatile pressure gradient values must therefore be increased. This would indicate that a greater pumping specification is needed in industrial applications. We note that for Fig. 3 the λ and Nf values imply a Darcy-Forchheimer flow. The influence of quadratic drag is quite low i.e. weak. Reynolds number is low so the Darcian drag is dominant in this regime and does exert a considerable influence on the velocity field.  Fig. 4 the influence of the porous matrix resistance embodied in the Darcian parameter λ (directly proportional to permeability) on the spatial distribution of velocity U versus Y is depicted. As a filtration material in various industrial pumping processes, λ is a key regulation parameter. Rising λ implies permeability increases and therefore the porous channel material in the channel e.g. soil, ceramic, insulation material, progressively decreases in presence i.e. less and less impeding material is present in the channel flow. We observe that as a result the Darcian bulk impedance is decreased substantially. The fluid receives less resistance and therefore is accelerated i.e. velocities increase as λ rises from 0.1 (low permeability i.e. densely-packed porous flter material) through 0.3, 0.5 to 1 and 2.
The profiles are symmetrically parabolic as in Fig. 2 since the channel prescribes no-slip conditions at either walls. Peak values of velocity arise in the centre line of the channel at Y = 0.5. Approximately a 20 % boost in U values accompanies a rise in λ from 0.1 to 2. Fig. 5 shows the velocity distribution versus Y for various Nf values i.e. Forchheimer quadratic parameter. A rise in N F from 1 (weak inertial effect) to 5, 10, 15 and 20, depresses the longitudinal velocity considerably. Quadratic inertial drag effects progressively decrease the influence of Darcian viscous-dominated drag. Velocity U is seen to fall from a peak value at the channel centre (Y = 0.5) for Nf = 1, to approximately 3.5 for Nf = 20 i.e. strong inertial drag. The porous filtration material in the channel has a significant deceleration effect on flow momentum at high Nf values. The influence of porous medium is infact much more effective than increasing the shearthickening i.e. rheological characteristic of the pulsating fluid (Fig. 3). Higher pulsation field values therefore serve to augment the filtration effect as the porous matrix secondorder drag is substantially depressed.
In Fig. 6 to 9, we provide the computations obtained for the special case of steady non-Newtonian flow in a Darcian channel. The influence of Reynolds number is illustrated in Fig. 6. The profiles correspond to β = 4, λ = 1 and Ps = 10. The velocity distributions are all skewed parabolas (to the right). The no-slip conditions (zero velocity) at each channel wall (Y = 0, 1) are observed as the zero U values at either end of the parabolas. Maximum velocities all occur at approximately Y = 0.65, unlike in the Newtonian case where they are exactly symmetrical. Velocity increases substantially with a rise in Reynolds number as this parameter rises from 0.5 (slow viscous flows) to 1, 2, 3 and 5. Up to Re ∼ 10 the Darcian model is valid and therefore λ is non-zero.
In Fig. 7 we have plotted the influence of rheological parameter (β) on longitudinal velocity. As in the general pulsatile, transient, rheological, non-Darcian case (Fig. 3) velocity U increases strongly with a rise in β. Hence the highest value of U (at the centre of the channel for each profile) increases from approximately 0.55 to 1.0 i.e. velocity is boosted by more than 80 %. We observe that in the absence of pulsating pressure gradient i.e. when no periodic pressure differences are present, the velocities in the channel are much lower than for the pulsating case (Fig. 3). As with Fig. 3 the fluid behaviour approaches Newtonian as β is increased; however the increase in β is much less effective than for the pulsating case (Fig. 3), verifying that in practical applications, pulsating pressure gradient is the dominant control parameter for achieving greater momentum transfer in the fluid.   8 illustrate the effects of the permeable material resistance (Darcian parameter λ) on the spatial distribution of longitudinal velocity U . As a characteristic of filtration materials in various industrial processes, λ is a key regulation parameter. As with the pulsating case (Fig. 4), we see that velocity rises markedly as λ rises from 0.1 (low permeability i.e. densely-packed porous filter material) through 0.3, 0.5 to 1 and 2. The fluid receives less resistance with rising Darcian parameter and therefore is accelerated. The profiles as in the pulsating case (Fig. 4) are once again symmetrically parabolic. Maximum values of velocity arise in the centre line of the channel at Y = 0.5. However velocities are much lower than in the pulsating case, once again illustrating the dominating influence of a periodic pressure gradient. Therefore we infer in consistency with other studies (e.g. Uchida [6], Wang [13] etc) that steady flows cannot attain anywhere near the velocity magnitudes achievable with pulsation mechanisms.
Finally in Fig. 9 we present the influence of steady pressure gradient component (Ps) on the flow velocity. Reynolds number has been fixed at unity. Velocity is seen to increase considerably as Ps rises from 10 through 12, 14, 16 and eventually to 20. Velocities are boosted from approximately 0.8 to 2.8 over this range. However we observe yet again that in the steady case, velocity magnitudes are much smaller than for any of the pulsatile flow cases (Figs. 3, 4 and 5). Therefore the presence of a pulsating pressure gradient is instrumental in generating high velocities whether the fluid is Newtonian or non-Newtonian. The present model is currently being extended to consider the numerical simulation of the pulsatile hydromechanics of more complex non-Newtonian fluids in porous media, including for example polar (couple stress) fluids and results will be communicated in the near future.

Conclusions
A finite element numerical solution for the pulsatile non-Newtonian flow in a porous medium channel has been presented. A non-Darcian drag force model has been employed to simulate the porous medium. The porous bulk matrix resistance (Darcian drag) is seen to substantially depress velocities i.e. velocities are increased for materials with progressively increasing permeability. Forchheimer effects also have a major decelerating influence on the flow velocity. Reynolds number is observed to largely boost velocity in the steady flow case, although profiles are seen to be skewed i.e. asymmetric. Velocities generally increase markedly with a decrease in non-Newtonian behaviour. The present computations also illustrate that much smaller velocities are generated in steady flow compared with pulsating flow.