Pulsating Flow of an Ostwald—De Waele Fluid between Parallel Plates

The flow between two parallel plates driven by a pulsatile pressure gradient was studied analytically with a second-order velocity expansion. The resulting velocity distribution was compared with a numerical solution of the momentum equation to validate the analytical solution, with excellent agreement between the two approaches. From the velocity distribution, the analytical computation of the discharge, wall shear stress, discharge, and dispersion enhancements were also computed. The influence on the solution of the dimensionless governing parameters and of the value of the rheological index was discussed.


Introduction
The laminar, oscillating flow, driven by a harmonic pressure gradient in Newtonian fluids in pipes has been studied theoretically at least since the work by Sexl in 1930 [1], who obtained the classical velocity profile in the radial direction in terms of the Bessel function of the first kind and order 0. In his paper, Sexl analyzed the behavior of the solution for some limit values of the dimensionless parameter n ⁄ , where is the radius of the tube, the angular frequency of the pressure gradient, and n the kinematic viscosity of the fluid. In 1952, Lambossy [2] re-analyzed the problem and explicitly concluded that the shape of the velocity profile depended on the dimensionless parameter n ⁄ . Later, Womersley [3,4] published a paper focusing on flows in arteries, in which it is stressed the importance of the dimensionless parameter used by Sexl and Lambossy, and that became known as Womersley number. The oscillatory flow of non-Newtonian fluids has been studied theoretically by Pipkin [5], Etter and Schowalter [6], and Rahaman and Ramkissoon [7] for viscoelastic fluids, among others. Uchida [8] studied the flow of a Newtonian fluid due to pulsatile pressure gradient (i.e., superposition of oscillatory and constant pressure gradients). Among the analytical studies on pulsatile flows of non-Newtonian fluids, the following works can be mentioned: Barnes et al. [9] for some kind of Oldroyd fluids; Phan-Thien [10,11] and Steller [12] for generalized Maxwell fluids; Davies et al. [13] for the Goddard-Miller model, powerlaw, and Segalman viscosity functions; Mai and Davis [14] for Bingham plastics; Daprà and Scarpi [15] for power-law fluids. Recognizing the non-Newtonian behavior of the blood [16], an important amount of work has been done addressing the circulation of blood in veins and arteries, stimulating the study of the pulsatile flow of non-Newtonian fluids (for example, references [17][18][19] among others). The most used approach to get the analytical solution of the non-Newtonian pulsatile flow is by means of power expansions, truncated at the first-or second-order, depending on the author. However, the increment of the flow rate with respect to the constant gradient pressure condition can be computed only with a second-order expansion [9][10][11]15]. Because most of the studies have been motivated by the transport of non-Newtonian fluids in the industry or by the blood motion in arteries and veins, the most used configuration corresponds to the flow in a cylinder, or with axial symmetry. Among the few analytical studies in non-cylindrical geometries, the works by Siginer [20] and Steller [12] can be mentioned. The first analyzed the pulsatile plane Pouseille flow of viscoelastic liquids of the memory integral type, and the second the flow of a generalized Maxwell fluid through a slit. Mc Ginty et al. [21] presented solutions of pulsatile flows of Newtonian, Maxwell, and Oldroyd B fluids in a cylinder and in the gap between two cylinders. Daprà and Sacarpi [22] analyzed the later configuration for a Bingham fluid. Letelier et al. [23] solved the pulsating flow of a linearly viscoelastic fluid in straight tubes of arbitrary cross-sections, presenting results of the flow for square and hexagonal ducts with rounded corners. Analytical solutions for two-dimensional configurations are rather scarce when compared with their cylindrical counterpart. In addition to the paper by Steller [12], we can mention those by Daprà and Scarpi [24], who got a perturbation solution to the secondorder of the pulsatile flow of a pseudoplastic Williamson fluid, and that by Nandakumar et al. [25] for a shear-thinning fluid in a two-dimensional stenosed channel.
Analytical solutions of pulsatile flows of non-Newtonian fluids are limited to very specific geometries, and they demand long algebraic developments. Thus, numerical approaches have been profusely used. The most used numerical approach to solve the equation of motions is the finite difference method. Among the articles that use this method for solving the pulsating and/or oscillatory flow of non-Newtonian fluids in cylindrical geometry, the following can be mentioned: Edwards et al. (1972, power-law fluid [26]); Balmer and Fiorina (1980, power-law [27]); Adusumilli and Hill (1984, truncated power-law [28]); Nakamura and Sawada (1987, a bi-viscosity model of a Bingham plastic [29]); Warsi (1994, power-law [30]); Mai and Davis (1996, Newtonian, power-law and a modified Bingham plastic [14]); Daprà and Scarpi (2006, power law [15]). Axisymmetric occluded pipes have been studied, for example, by Javadzadegan et al. (2009, Oldroy-B and Cross model [31]). Rectangular or two-dimensional geometries have been studied by Yakhot et al. (1999, Newtonian fluid [32]) and Daprà and Scarpi (2007, Williamson fluid [24]). Flow in the gap of coaxial cylinders of a Bingham fluid under periodic pressure gradient and periodic oscillation of the inner cylinder has been studied by Daprà and Scarpi [22]. The work by Tu and Deville [33] is among the exceptions in which the finite element method is used to solve the pulsatile flow of power-law and Herschel-Bulkley fluids in a cylindrical pipe with an occlusion. Pontrelli [34] complemented an implicit finite difference method in time with a spectral collocation method in space in order to solve the pulsatile flows of fluids modeled as Oldroy-B and an ad hoc blood model.
Most of the experimental studies concerning pulsating flows of non-Newtonian fluids are oriented to determine the enhancement of the temporal mean discharge when compared with the constant pressure gradient driven flow. Some of them also deal with the effect that pressure pulsation has on energy consumption. The motivation of most of the papers comes from studies on blood circulation. Barnes et al. [9] got experimental relationships between the flow enhancement, the frequency of the fluctuation, pressure gradient, and amplitude of the oscillation for Newtonian and shear-thinning fluids generated by an aqueous solution of polyacrylamide. Harris and Maheshwari [35], also with solutions of polyacrylamide, compared the predicted and measured velocity profiles. Their work did not use specific rheological models, but they were within the general class up to second-order viscoelastic fluids.
Thurston [36], in 1975, performed measurements of pressure and volume flow of blood considering steady, oscillatory, and pulsatile flows. The blood is considered a viscoelastic fluid, and the comparison with measurements on a Newtonian glycerol solution helped the authors to specify the special features of the blood flow. Davies and Chakrabarti [37] worked with solutions of polyacrylamide in order to experimentally obtain the enhanced discharge and energy consumption and compared them with the corresponding values for Newtonian fluids. Along the same line, Phan-Thien and Dudek [10] from their experiments with aqueous solutions of Separan AP-30 found that the enhancement decreased with the frequency. Kajiuchi and Saito [38] performed experiments with clay slurries that behaved as Bingham plastics, finding that the flow enhancement could be expressed as a function of the Stokes number for a Bingham plastic fluid, the Bingham number, and the ratio between the pressure oscillation amplitude and the gradient of the base pressure. Lin et al. [39] studied non-colloidal suspensions modeled as Newtonian fluids for a shear rate below a cut value and as power-law fluids for shear rates greater than that value. Among other conclusions, they determined that particle migration might significantly affect flow enhancement.
Most of the researches involving pulsating flows deal with cylindrical or distorted cylindrical geometries because they have been motivated by studies of blood flows in veins and arteries. However, the flow in fractured rocks entails fluid motion through slender gaps that can properly be considered as the flow between parallel plates. The pulsatile flow of shear-thinning fluids through fractured rocks is found in enhanced oil recovery operations [24].
The goal of the present study was to determine analytical solutions for the two-dimensional flow of an Ostwald-de Waele type fluid driven by a pulsating pressure gradient. To accomplish this goal, a perturbation method was used to obtain the velocity distribution up to a second-order term, from which the instantaneous discharge, wall shear stress, cycle-average discharge, and dispersion coefficient were also analytically obtained and computed. A numerical solution using a spectral method was computed in order to check the analytical solution of the velocity distribution. As far as the authors are aware, no solution for the two-dimensional pulsatile flow of an Ostwald-de Waele fluid has been previously published. Daprà and Scarpi [15,24] solved this problem in cylindrical coordinates and used a two-dimensional geometry for a Williamson fluid. The solution presented in this article includes more details of the algebra than reported previously [15,24].

Analytical Solution
The laminar, incompressible, two-dimensional, pulsatile flow between two parallel plates of a power-law fluid is governed by the continuity and momentum equations, which are reduced to: where is the pressure, is the velocity component in the direction, parallel to the plates. The axis is normal to the plates, with the origin located at the inferior plate. is the time. , , and are the density, consistency coefficient, and flow behaviour index of the fluid, respectively. In pulsatile flow, the pressure gradient is given by: where and are the frequency and the relative amplitude of the oscillating part of the pressure gradient, respectively. For infinitely long plates, the velocity is independent of the variable . Two parallel plates separated at a distance 2 have a symmetry plane along = and non-slip conditions at the plates. Thus, the boundary conditions are Assuming > 0 and ≪ 1 , the velocity derivative satisfies ⁄ ≥ 0 for 0 ≥ ≥ , and Equation (2) becomes: Denoting with an upper asterisk the dimensionless variables: * = , * = , * = , The velocity scale is not arbitrary and has been defined later. Dropping the asterisks of the dimensionless variables, the momentum equation is written as: where W and l are dimensionless parameters defined as: As there is not a velocity scale imposed externally to the problem, should depend on the other parameters, i.e., , , and , resulting = ⁄ . Thus, the dimensionless frequency becomes Replacing in terms of , it is easy to recognize that l corresponds to the inverse of a Reynolds number associated with power-law fluids, l = 1 ⁄ = ⁄ . In terms of : Equation (8) becomes: where n denotes the apparent viscosity: n = , with the boundary conditions: The above differential equation will be solved expanding it in a series of up to the secondorder: Replacing the expansion given by Equation (16) in Equation (13) yields to: The solution at the zero-order is the classical result for a power-law fluid, and it is obtained after integrating the following equation with the boundary conditions: resulting the zero order velocity distribution: The unsteady characteristic of the flow is considered from the first-order onwards. Using the solution for , the differential equation for order is Considering that the first-order solution of the velocity has the form: where Â [ ] is the real part of the complex , and replacing in Equation (22), the differential equation becomes Expanding the derivative and rearranging terms, the above equation is written as: Equation (25) is a non-homogeneous differential equation. The function is split into two components, = + , where is the solution of the homogeneous part, and = −1 W ⁄ the solution of the non-homogeneous one. The differential equation for is given by: where the new variable h is defined as: Equation (26) is an Emden-Fowler type differential equation, with solution where (n) is the Bessel function of the first kind and order n. and are the integration constants. and n are defined as Imposing the boundary conditions (Equations (19) and (20)), the values of and are obtained: Thus, the first-order term of the velocity is: The second-order involves the solutions of zero and first-order. Using the complex identity The velocity is assumed to be composed by a steady and an unsteady component: Thus, the steady-state component is obtained from Using | | = , ( ) = ( ̅ ), the expression for the series expansion of the product of Bessel functions (reference [40], page 148), the term , can be written as Thus, the differential equation for is given by where = ( − 1) 4 Given the structure of Equation (38), a series solution for is feasible, resulting: where the coefficients are given by: Obviously, the solution given by Equation (40) satisfies the boundary conditions. The differential equation for the unsteady part of the second-order solution is: The structure of is similar to that for the unsteady part of the first-order solution, i.e., is obtained from: After a rather cumbersome algebra: Finally, the second-order term of the velocity is: Thus, inserting the solutions given by Equations (21), (31), and (49) into Equation (16), the velocity distribution up to the second-order is obtained. The shear stress is made dimensionless with , and it is computed directly from = 1 (50) After some algebra, the second-order approximation of is found to be: from where the shear stress acting on the bottom ( = 0) is determined: It is easy to show that the bottom shear stress averaged on a cycle is equal to the value associated with the steady flow, i.e., 〈 〉 = 1, where the angular bracket denotes cycle average.

Flow Enhancement
The instantaneous dimensionless discharge can be expressed as The flow enhancement is defined as = (〈 〉 − ) ⁄ , where 〈 〉 is the average discharge of the pulsatile flow along a period of fluctuation, and is the discharge associated with the steady component of the pressure gradient. In this problem, = . Contribution of is nil, resulting in that the flow enhancement is only due to . The second-order discharge (Equation (56)) has two addends, but as the second one depends on , it does not contribute to the increment of discharge in a period. Thus, the flow enhancement is given by: In the above equation, is obtained from Equation (41). The analytical expression for the flow enhancement when inertia is negligible can be obtained, making zero the left-hand-side of Equation (8) and integrating along , resulting = . This is the same discharge increment than for a pulsating flow in a cylindrical pipe [37]. The dimensionless cycle-averaged cross-section-mean-velocity is easily computed from Equation (53), resulting 〈 〉 = (2 (2 + 1) ⁄ ) (1 + (1 − ) (4 ) ⁄ ).

Dispersion Coefficient
The dimensionless coefficient of longitudinal dispersion, * , is defined as ( [41], p.86 where is the dispersion coefficient, is the coefficient of molecular diffusion, and is the dimensionless velocity deviation from the vertical average ( ). The structure of the triple integral given by Equation (58) precludes us to look for a second-order analytical solution for * , but it is possible to compute its cycle-average value. The dimensionless dispersion coefficient averaged in a period of oscillation is 〈 * 〉 = ∫ * d . The order of integration can be changed using the Fubini theorem [42], and the cycle average is computed from The velocity deviation contains terms involving and , whose integrals in a period are nil. Thus, it is not necessary to carry out the triple integrals in for many of the terms that appear in Equation (58). Then, the inner integral is: From the above result, it is easy to see that the cycle averaged dispersion coefficient has the structure 〈 * 〉 = 〈 * 〉 + 〈 * 〉 (60) with 〈 * 〉 = 2 3(2 + 1) (4 + 1)(5 + 2) and After some rather boring algebra, the term accompanying can be computed, and the cycle averaged dispersion coefficient is found to be: where , = (2 + 1)(2 + 2 + 8 + 2) (2 + 2 + 2 + 1)(2 + 2 + 4 + 1)(2 The dispersion enhancement is given by * = (〈 * 〉 − 〈 * 〉 ) 〈 * 〉 ⁄ . Because the intricacy of the coefficient of second-order, it is not possible to get a simpler expression like that obtained for the discharge (Equation (57)).

Numerical Solution of the Velocity Distribution
In order to verify the correctness of the analytical expression of the velocity profiles reported in the previous section, Equation (13) has been numerically integrated with a spectral method. Spectral methods constitute the most confident approach to the solution of fluid mechanics equations for both Newtonian (see for instance [43][44][45]) and non-Newtonian (see for instance [46][47][48]) fluids, reducing at the minimum the impact of the numerical scheme on the physics of the flow under investigation. Moreover, the higher computational cost of these methods is widely compensated by their fast (exponential) rate of convergence [49]. The spatial discretization of Equation (13), along with the boundary conditions (14)- (15), has been carried out through a spectral patching collocation method based on Gauss-Lobatto-Chebyshev quadrature nodes (see for instance [50]). Denoted with , with = 1, . . where is the -term of the Chebyshev collocation derivative matrix [49], and h the value of the apparent viscosity pertaining to the grid point. In order to deal with the unbounded shear stress gradients associated with the non-Newtonian power-law model, we have adjusted the power-law model by replacing the apparent viscosity n as follows [51][52][53]: with an adjustment parameter for which non-negative values have to be prescribed. When = 0, the standard power-law model is recovered. The influence of the value on the results has been discussed later on. The homogeneous Dirichlet and Neumann boundary conditions (14) and (15) Figure 1 shows that the results are essentially independent of the value. Similar results (not shown herein) have been obtained, varying each parameter of the ( , W) pair in the range (1,100). In what follows the value = 10 has been considered. Figure 2 compares the velocity profiles analytically deduced with the corresponding ones numerically calculated. Four instants have been considered, namely, = 0.25, 0.5, 0.75, 1.0 . The power-law exponent is equal to 0.5, while the Reynolds number, the pulsation Ω , and the coefficient are 1.0, 1.0, and 0.1, respectively. Figure 2 shows that the analytical solution is in very good agreement with the numerical one, for all instants.  = . = .
As expected, similar conclusions cannot be drawn when the oscillation amplitude is increased.
Indeed, Table 3, where both the percentual errors in the case = 0.5 = 1, Ω = 1 , and = 0.5 are reported, shows that as the perturbation parameter is increased, the accuracy of the proposed solution, as it could be expected, significantly reduces.

Wall Shear Stress
The wall shear stress is given by Equation (52). Due to the cumbersome structure of the relationship, it was not feasible to determine analytically some characteristics like the mean, maximum, and minimum values of the wall shear stress. Some limit behaviours of were analyzed for values of the flow index = 0.5, 0.75, and 1 (shown below).

Wall Shear Stress for W ≫ 1 and ~1
For large dimensionless frequency and ~1, Equation (13) indicates that W ⁄~ ⁄ , i.e., the contribution of the fluctuating part was negligible, as it is corroborated in Figure 3, where the values W = 100, = 1, and = 0,1 were used. It was observed that the amplitude deviated at the most 1% from the non-pulsating condition (for = 1). Also, it showed a delay of the response with respect to the pressure fluctuation. The non-linearity reduced even more the pulsation, already reduced by the effect of large W, with a decrease in amplitude as the flow index decreased.

Wall Shear Stress for W ≪ 1 and ~1
For small W and ~1 , the terms of first-and second-order in the solution of became important as they were divided by W, resulting in higher amplitude fluctuations, as it is observed in Figure 4, computed for W = 0,1, = 1, and = 0.1.

≪ 1 and W ~1
For small and W ~1, the term of order in Equation (52) could not be neglected, and the pulsation became important, with little influence of the flow index, as could be seen in Figure 6, where the three curves were practically coincident for W = 1, = 0,1 and = 0.1.

Discharge
The discharge was computed from Equations (52) The above result indicated that the dependency of the velocity profile with time was negligible, from where no change of the discharge with time was expected. As an example, Figure 7

Discharge for W ≪ 1 and ~1
When W was small, the pulsating term was not negligible, and the fluctuation in time was noted in the discharge. Considering W ~ ≪ 1 resulted that the discharge was the superposition of the discharge for the steady-state condition with a fluctuating behaviour of the velocity given by , resulting in a harmonic oscillation of the discharge around the steady flow. As an example, Figure 8 shows the discharge for W = 0.1, = 1, = 0.1 for = 0.5, 0.75, and 1.  When W ~1 and ≪ 1, the contributions of the unsteady components of the discharge became important and could not be neglected as in the previous case. As an example, Figure 10

Flow Enhancement
The increment of discharge due to the non-linear behaviour of the constitutive law, summarized by (Equation (57)), is exemplified in Figures 11 and 12. The first figure presents the flow enhancement for = 0.75 as function of the dimensionless frequency with the Reynolds number as a parameter. In Figure 12, for = 1, the flow index was considered as a parameter. In both cases, = 0.1. As expected, at low frequencies, the discharge increment went to the value determined analytically. At large frequencies, there was no enhancement with respect to the non-pulsatile condition. Although the value of W, at which was nearly zero, increased with , the product W remained nearly constant and was equal to 100. The same happened with the value of W at which the enhancement reached its maximum value. It decreased with , but W ≈ 0.1. For a given , the W span, at which decreased from its maximum value to zero, was about 3 decades. As expected, the closer the fluid to a Newtonian behavior, the lower the flow enhancement, ( → 1, → 0), as shown in Figure 12.

Dispersion Enhancement
The dispersion enhancement, characterized by , behaved in an identical manner as , as shown in Figures 13 and 14, where the plots are for the same values of and compared to in Figures 11 and 12. Thus, it was observed that the dispersion enhancement goes to zero at a large dimensionless frequency as the pulsating pressure effect become negligible. For small frequencies, tend to a constant value. Like , the range of W, at which decreases from its maximum value to zero, was around 3 decades, satisfying W ≈ 0.1 for the maximum and W ≈ 100 for zero, with → 0 as → 1 (Figure 14).

Conclusions
The pulsating two-dimensional flow of an Ostwald-de Waele fluid was solved analytically using the perturbation method up to the second-order term. Details of the analysis were presented, skipping only the more cumbersome algebra. A numerical solution based on spectral methods was also used to compute the velocity profile. Both analytical and numerical solutions showed very good agreement. From the velocity distribution, analytical expressions for the discharge and the wall shear stress were determined. The discharge and dispersion enhancements resulting from the nonlinearity of the rheology were also determined. These enhancements were noticed at the second-order. In terms of the dimensionless frequency, the enhancement coefficients behaved in a similar way for both the discharge and the dispersion coefficients, presenting a variation along 3 decades of W from their maximum values (corresponding to the non-inertial limit) to zero (for the non-pulsatile condition).