Shear-Thinning Fluid Flow in Variable-Aperture Channels

: Non-Newtonian fluid flow in a single fracture is a 3-D nonlinear phenomenon that is often averaged across the fracture aperture and described as 2-D. To capture the key interactions between fluid rheology and spatial heterogeneity, we adopt a simplified geometric model to describe the aperture variability, consisting of adjacent one-dimensional channels with constant aperture, each drawn from an assigned aperture distribution. The flow rate is then derived under the lubrication approximation for the two limiting cases of an external pressure gradient that is parallel/perpendicular to the channels; these two arrangements provide upper and lower bounds to the fracture conductance. The fluid rheology is described by the Prandtl–Eyring shear-thinning model. Novel closed-form results for the flow rate and hydraulic aperture are derived and discussed; different combinations of the parameters that describe the fluid rheology and the variability of the aperture field are considered. The flow rate values are very sensitive to the applied pressure gradient and to the shape of the distribution; in particular, more skewed distribution entails larger values of a dimensionless flow rate. Results for practical applications are compared with those valid for a power-law fluid and show the effects on the fracture flow rate of a shear stress plateau.


Introduction
Non-Newtonian fluid flow in fractured media is of interest for many environmentally related applications, such as hydraulic fracturing, drilling operations, enhanced oil recovery, and subsurface contamination and remediation.The basic building block in fractured media modeling is a thorough understanding of flow and transport in a single fracture [1].A key concept in single-fracture flow and transport is the fracture aperture, defined by the space between the fracture walls.As a result of the heterogeneity of these surfaces, the fracture aperture is spatially variable.To model this variability, two basic approaches have been adopted.The first describes the aperture as a 2-D random field, typically described by an autocorrelation function of finite variance and integral scale [2].The second envisages the aperture variability as the outcome of the joint variation in self-affine surfaces, correlated at all scales [3].In both cases, flow modeling at the single-fracture scale leads to the determination of the flow rate under a given pressure gradient as a function of the parameters that describe the variability of the aperture field or of the confining walls.A hydraulic aperture can then be derived from the flow rate [4] as the aperture of a smooth-walled conduit that would produce the same flow rate under a given pressure gradient as the real rough-walled fracture.
When fluid behavior is non-Newtonian, the effects of spatial variability are compounded with the influence of rheology, producing striking results, such as pronounced channeling effects [5].Different constitutive equations have been used to represent non-Newtonian behavior in fracture flow, ranging from the simpler two-parameter power-law model [5] to the four-parameter Carreau-Yasuda equation [6].A comprehensive comparison of results for different constitutive equations is still lacking, but the impact of fluid rheology is likely to be significant in view of the diverse behavior of rheological equations in the vicinity of the zero-shear-rate limit [7].
Detailed 2-D or 3-D flow modeling of non-Newtonian flow in single fractures needs to be tackled numerically, with considerable computational effort, given the nonlinearity of the flow.Not surprisingly, some authors have pursued a simpler approach, with the aim of providing order of magnitude estimates and reference benchmarks for fracture conductivity.Basically, this approach considers a simplified, extremely anisotropic fracture geometry, with the aperture variable along one direction and constant-aperture channels along the other.The arbitrary orientation of the external pressure gradient with respect to the channels gives rise to two limit cases: (i) the parallel arrangement, which provides an upper bound to the conductivity, and (ii) the serial arrangement, which provides a lower bound.Flow in an isotropic aperture field is then addressed by considering the fracture as a random mixture of elements in which the fluid flows either transversal or parallel to the aperture variation.The hydraulic aperture is derived by a suitable averaging procedure, previously adopted for a deterministic aperture variation [8].
The present paper follows this avenue of research by exploring the impact of a classical, twoparameter shear-thinning constitutive equation, the Prandtl-Eyring model [9], which overcomes the unrealistic behavior of the power-law model, having infinite apparent viscosity for the zero shear rate.
Several suspensions of practical interest, such as magnetohydrodynamic (MHD) nanofluids, are well-interpreted by the Prandtl-Eyring model [10].Many nanofluids are adopted in technologies such as fuel cells, microelectronics, and hybrid-powered engines.Nanofluids are used in industrial technologies because of their utility in the production of high-quality lubricants and oil, as they flow in fractures and small channels.
Section 2 derives the flow rate under an assigned external pressure gradient for the flow of a Prandtl-Eyring fluid in a parallel-plate fracture.Section 3 presents the simplified geometry adopted, derives general expressions of the flow rate for flow that is parallel or perpendicular to constantaperture channels, and proposes a method to evaluate the hydraulic aperture for the 2-D case.Section 4 introduces a specific probability distribution function for the aperture-the gamma distributionand illustrates results for the flow rate and hydraulic aperture for the two different 1-D cases.A dimensional comparison with results for the power-law fluid is drawn in Section 5. Conclusions and suggestions for future work are listed in Section 6.

Prandtl-Eyring Fluid Flow in Constant-Aperture Fracture
We consider the flow of a non-Newtonian Prandtl-Eyring fluid between two smooth parallel plates separated by a distance b (fracture aperture); the coordinate system is shown in Figure 1, with upper and lower plates at  = ±  2 ⁄ .A uniform pressure gradient  = −  ⁄ is applied in the direction, where the generalized pressure  =  +  includes gravity effects,  is pressure,  is gravity, and  is fluid density.Assuming flow in the -direction, the velocity  is solely a function of  under the hypothesis of infinite width W. Momentum balance yields a linear shear stress profile: A Prandtl-Eyring fluid is described rheologically in simple shear flow by [9]: where the parameter  = 1  ⁄ of dimensions (T −1 ) is the inverse of relaxation time  , and the parameter  is the Eyring characteristic shear tress of dimension (ML −1 T −2 ).For a vanishing shear rate (  ⁄  → 0), the behavior tends to Newtonian with viscosity  = /, as is easily seen via the first-order expansion sinh  ≅ .The behavior for a nonzero shear rate is shear-thinning, with a vanishing apparent viscosity for high shear stress.Figure 2 shows the apparent viscosity ( ), defined by the relationship  = ( ) , as a function of the shear rate  , for realistic parameter values.Substituting Equation (2) in Equation ( 1) and integrating with the no-slip condition at the wall  (∓  2 ⁄ ) = 0 gives the velocity profile between  = −  2 ⁄ and  = +  2 ⁄ as The total flow rate  through the fracture for a width  in the y-direction perpendicular to the pressure gradient is derived by integrating Equation (2); the result is where  is the flow rate per unit width, and  is the average velocity.For the vanishing dimensionless pressure gradient (  (2) ⁄ → 0), the hyperbolic functions can be replaced with their second-order expansions cosh  ≅ 1 +  2 ⁄ and sinh  ≅  +  6 ⁄ , and the flow rate becomes which is the "cubic law" [11] valid for a Newtonian fluid.
If the aperture varies, as in real rock fractures, a flow law of the type in Equations ( 4) and ( 5) is valid in replacing the constant aperture  with a hydraulic (or flow or equivalent) aperture  , accounting for the variation in fracture aperture [8,12,13].The availability of a closed-form expression for the hydraulic aperture as a function of fluid rheology and aperture variability is of particular interest in view of its use in numerical simulators.

Flow in Variable-Aperture Channels
Flow and transport simulations in variable-aperture fractures typically consider either the aperture (, ) to vary as a two-dimensional, spatially homogeneous, and correlated random field with a probability density function () and assigned statistics, or to have walls described by a fractal distribution of a given Hurst coefficient H, correlated at all scales [3].In the former case, the fracture dimensions are assumed to be much larger than the integral scale of the aperture autocovariance function; then, under ergodicity, spatial averages and ensemble averages are interchangeable, and a single realization can be examined [4].This approach was followed by Di Federico [14] and Felisa et al. [15] to study the flow of power-law and truncated power-law fluids in simplified aperture fields in which the aperture varies only along one spatial coordinate, and the external pressure gradient (and hence the flow) is either transverse or parallel to the aperture variability; such an idealized fracture of dimensions  and  is shown in Figure 3. Here, we consider the flow of a Prandtl-Eyring fluid in both arrangements.

Channels in Parallel
Consider flow along the direction  parallel to constant-aperture channels and driven by the external pressure gradient  , as illustrated in Figure 3a.To obtain the volumetric flux, a procedure similar to that adopted by Zimmerman and Bodvarsson [16] is used.The fracture, which has length  and width , is discretized into N neighboring parallel channels, which are all of equal width  =   ⁄ , length  , and constant aperture  .Assuming that the shear between neighboring channels and the drag against the connecting walls may be neglected, the flow rate in each channel of a constant aperture along the -direction is given by Equation ( 4) with  in place of  and  in place of .Hence, summing over all channels, the total flow rate in the fracture is Taking the limit as  → ∞, the width of each channel tends to zero, and the discrete aperture variation tends to a continuous one; then, under ergodicity, Equation (8) gives the following expression for the flow rate per unit width: where () is the probability distribution function of the aperture field , defined between 0 and ∞.Finally, the hydraulic aperture  may be derived numerically by equating  from Equation (9) to Equations ( 4) and ( 5) written with  in place of .
The assumption of negligible shear between neighboring channels in the limit As a first approximation, the shear rate in the x-y plane due to a varying b is with a shear stress (see Equation ( 2)) equal to Substituting into Equation (10) yields For smooth variation in b, with 0 and 0, the condition in Equation ( 13) is satisfied since, on the right-hand side, we have a finite-order term that is always >4, and on the left-hand side, we have a small-order term.

Channels in Series
Consider flow along the direction  that is perpendicular to constant-aperture channels and driven by the external pressure gradient  , case (b).The fracture, which has length  and width , is discretized into N cells in series of equal length  =   ⁄ , each of width  and constant aperture  .As volumetric flux  through each cell is the same, the flow rate per unit width in each cell is also the same, i.e.,  =  =   ⁄ .The total pressure loss along the fracture in the -direction, ∆, can be expressed as the sum of pressure losses in each cell, ∆ , as having neglected the pressure losses due to the succession of constrictions and enlargements.In turn, dividing by  yields the external mean pressure gradient  as where the pressure gradient in each cell of constant aperture  is given by  =   ,  , ,  , which is obtained by deriving the pressure gradient as a function of flow rate from Equations ( 4) and ( 5), written by replacing the subscript  with  and  with  .Taking the limit as  → ∞, the length of each cell tends to zero, and the discrete aperture variation tends to a continuous one; then, under ergodicity, Equation (15) gives the following expression for the mean pressure gradient in the -direction: The integration of Equation ( 16) implicitly gives the flow rate as  =   , , , () .Finally, the hydraulic aperture  may be derived numerically upon equating  derived with Equation ( 5) written with the subscript  in place of  and with  in place of .
As for the parallel case, in the serial arrangement, it is assumed that which is satisfied for smooth variations in b.Abrupt variations in b could be modeled by inserting a local dissipation due to, e.g., an expansion of the flow.

Flow in 2-D Isotropic Aperture Field
Non-Newtonian flow in a fracture that is characterized by an isotropic, two-dimensional aperture variation is highly complex, as shown by, e.g., De Castro and Radilla [6], and the hydraulic aperture can be obtained only by means of numerical simulations.However, it can be argued [16] that the scheme with channels in parallel is an upper bound to the hydraulic aperture for the general 2-D case, while the scheme with channels in series provides a lower bound, in analogy with the wellknown expressions for hydraulic conductivity [17].If the fracture is seen as a random mixture of elements in which the fluid flows either transversal or parallel to aperture variation, the flow can be approximated by a suitable average of these flows; the ergodicity assumption ensures that boundary effects are negligible [4].Hydraulic aperture values derived for the two schemes differ significantly; hence, following the procedure adopted in the literature [4,8,12,14,15], an estimate of the hydraulic aperture is derived as the geometric mean of the hydraulic apertures for the parallel and serial arrangement as  =   . (18)

Aperture Probability Distribution
A gamma distribution of the shape parameter  and scale parameter  , entailing non-negative apertures, is adopted to quantify Equations ( 9) and ( 16), consistent with earlier work [15].Its probability density function, expected value, variance, and skewness are given by where Γ() is the gamma function.The gamma distribution is illustrated in Figure 4; for  = 1, it reduces to the exponential distribution with maximum skewness, while as  → ∞, the gamma distribution tends to a normal distribution with the same mean and variance and zero skewness.

Channels in Parallel
Inserting Equation ( 19) into Equation ( 9) gives the following result after integration [18] (p.403), some algebraic manipulations, the use of Equation (20), and the exploitation of the properties of the gamma function: where The relationships in Equations ( 24)-( 26) establish that the ratio between the shear stress  at the wall of a parallel-plate fracture of aperture 〈〉 and the characteristic shear stress  describing the Prandtl-Eyring fluid cannot exceed the shape parameter  of the aperture distribution.Figure 5 shows the dimensionless flow rate per unit width  as a function of Ω for different values of .It is seen that as Ω → 0 , the flow rate tends to −∞; such negative values are obviously not realistic.On the other hand, for Ω → , the flow rate tends to infinity, and the curves show a vertical asymptote of the equation Ω = .Similarly, such results are unrealistic.For low values of Ω, the curves tend to overlap, regardless of the value of .In general, the flow rate depends, in a nonlinear fashion, on the dimensionless pressure gradient and the distribution parameters.The hydraulic aperture for the parallel arrangement  is obtained by solving the following implicit equation in the unknown  : Figure 6 illustrates the ratio  of the hydraulic to mean aperture defined by Equations ( 27) as a function of Ω for different values of .The ratio  strongly increases with Ω, more so for lower values of  , i.e., a more skewed distribution.The curves show a vertical asymptote when the dimensionless pressure gradient Ω approaches  .For low values of Ω , the curves are almost horizontal.

Channels in Series
The pressure gradient cannot be obtained explicitly as a function of flow rate from Equations ( 4) and ( 5) written by replacing the subscript  with  ; hence, the integral in Equation ( 16) with Equations ( 19)-( 21) can only be solved numerically.Doing so yields  =   , , , 〈〉,  , which can be made dimensionless as  =  (Ω, ), where  is defined by Equation ( 23) with the subscript  in place of , and Ω is defined by Equation (24).
The hydraulic aperture  is obtained by solving the following implicit equation in the unknown  :

Results for a Power-Law Fluid
A power-law fluid has the rheological equation where  is the rheological index, and  is the consistency index of dimensions (ML −1 T n−2 ); for  < 1 and  > 1 , the model describes shear-thinning (pseudoplastic) and shear-thickening (dilatant) behavior, whereas for  = 1, Newtonian behavior is recovered, and  reduces to dynamic viscosity .The relationship between the applied pressure gradient  and flow rate per unit width  in a parallel-plate fracture of aperture  is For the gamma distribution of an aperture defined by Equations ( 19)-( 22), results are obtained after some algebra as which further reduces for a Newtonian fluid to  =  = 1 12 ⁄ .

Comparison between Power-Law and Prandtl-Eyring Models
In this section, the parallel arrangement is considered, and a comparison between the powerlaw and the Prandtl-Eyring models is proposed.Here, we aim to analyze the effects of a shear stress plateau in the Prandtl-Eyring model for low shear rates.The nature of a non-Newtonian fluid is hard to model, and this inevitably leads to a complex expression for the fluid rheology.The most commonly used model for a shear-thinning fluid is the power law relationship; it is a valid tool for studying non-Newtonian fluid behavior without introducing excessive mathematical complexities.It is, however, not entirely realistic, as for low shear rates ( → 0), the apparent viscosity tends to infinity ( → ∞), while for ( → ∞), it results in  → 0. The Prandtl-Eyring model is characterized by a viscosity upper bound; this rheologic model presents a plateau, and for  → 0, the fluid resembles a Newtonian fluid with  → /.To compare the two rheologies, a Prandtl-Eyring fluid ( = 0.003465 Pa and  = 0.0231 s ) is considered [7] and adopted to evaluate the parameters of a comparable power-law fluid by means of a least-squares fit performed in the range of a shear rate between 10 and 1 × 10 3 s −1 : the result of the best fit yields  = 0.0175 Pa × s and  = 0.13 (Figure 7).To analyze the influence of the plateau, we define the shear rate ratio 〈 〉 /〈 〉 , where 〈 〉 is the average fracture shear rate of a Newtonian fluid of dynamic viscosity  = /, written as and similarly, the average fracture shear rate of a Prandtl-Eyring fluid is Figure 8 depicts the flow rate versus pressure gradient for power-law and Prandtl-Eyring fluids for channels in parallel within the range of validity (Equation ( 26)) of Equation ( 23).Both the red and blue solid lines, representing the Prandtl-Eyring and the power-law fluids, respectively, are monotonically increasing functions of the pressure gradient, with the former model showing a flow rate several orders of magnitude higher than the latter.The power-law model presents a linear trend in log-log coordinates, while two different behaviors can be observed in the Prandtl-Eyring model: for low-pressure gradients, the slope of the red solid line is minor compared with the one of the blue solid line but higher than the Newtonian case; on the other hand, when the average fracture shear rates of the two fluids are sufficiently different, the Prandtl-Eyring model rapidly increases its steepness.When the average fracture shear rates 〈 〉 tend to 〈 〉 , most of the channels undergo a flow regime of low shear rates, where the Prandtl-Eyring model simply resembles a Newtonian fluid of dynamic viscosity  = /.

Conclusions
Non-Newtonian flow in a fracture characterized by an isotropic, two-dimensional aperture variation is complex; a simplified view is proposed that envisages the fracture as a random mixture of elements in which the fluid flows either transversal or parallel to a 1-D aperture variation, which is described by a probability distribution function of assigned mean and variance.The Prandtl-Eyring model was chosen to describe the fluid, as it captures the shear-thinning behavior of most fluids of interest in fractured media flow without the unrealistically large apparent viscosity for low shear rates typical of power-law fluids.The latter model was chosen for comparison.The gamma distribution was selected to represent aperture variability.Key results are as follows:

•
Values of the flow rate are extremely sensitive to the applied pressure gradient and to the shape of the distribution; in particular, a more skewed distribution entails larger values of the dimensionless flow rate; • A comparison was drawn between the Prandtl-Eyring (PE) and the power-law (PL) model having equal apparent viscosity for a wide range of shear rates.For channels in parallel, the absence of a shear stress plateau for low shear rates, associated with power-law rheology, implies an underestimation of the fracture flow rate with respect to the Prandtl-Eyring case.For the latter fluid, low-pressure gradients are characterized by a flow regime dominated by the plateau viscosity  = / , with a behavior similar to a Newtonian (N) fluid, while for sufficiently high-pressure gradients, as the ratio between the two average fracture flow rates (PE

Figure 4 .
Figure 4. Gamma distribution pdf for different values of the parameter  and assuming a mean fracture aperture 〈〉 = 1 mm.

Figure 5 .
Figure 5. Dimensionless flow rate per unit width  versus dimensionless pressure gradient Ω for different values of the distribution parameter .

Figure 6 .
Figure 6.Ratio  =  〈〉 ⁄ versus dimensionless pressure gradient Ω for different values of the shape parameter .

Figure 8 .
Figure 8. Plot of flow rate vs. pressure gradient; the red solid line represents the Prandtl-Eyring model ( = 0.003465 Pa and  = 0.0231 s ), the blue solid line represents the power-law model ( = 0.0175 Pa × s and  = 0.13), and the black dashed line represents the Newtonian model ( = / = 0.15 Pa × ).The plot of shear rate ratio 〈 〉 /〈 〉 vs. pressure gradient is illustrated by the black solid line.