Original Research ARTICLE
Effective rheology of Bingham fluids in a rough channel
- 1CNRS, FAST, Université de Paris-Sud, Orsay, France
- 2Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway
We derive the volumetric flow rate vs. pressure drop of a Bingham fluid in one-dimensional channels of variable apertures in the lubrication approximation. A characteristic length scale, a* characterizing the flow is introduced in order to distinguish between a high and a low flow rate regime. We illustrate the calculation for channels with periodically varying apertures. We then go on to consider apertures that are self affine. We determine how the scaling properties of the aperture field is reflected in the effective flow equations. Finally, a series expansion for high and low flow rates that works very well over the entire range of flow rates is proposed. This truncated expansion allows us to predict the domain of validity of the two expansions by comparing a* to the aperture distribution.
Yield-stress fluids, viz. fluids that require that the applied stress is above a non-zero threshold—a yield stress—for flowing , are found in many practical applications, such as in oil industry. Indeed, emulsions , mud [3, 4], heavy oil , polymeric gels such as carbopol  and foams  generate non-zero yield stress for flowing. In an oil recovery context, such fluids are injected to control premature production of water or gas. They do this by reducing poor sweep efficiency due to reservoir heterogeneities [7, 8]. Hydraulic fracturing operations—now very much a timely topic because of the increasing importance of fracking—involve also the flow of cross-linked polymer gels that are characterized by a yield stress to help carry proppant particles through the fractures. Common for these examples is that they occur in complex geometries; in the first in porous media, whereas in the second, the flow occurs inside rough cracks.
There have been a number of experimental [2, 9], numerical [10–16] and theoretical [17, 18] studies of the flow of yield stress fluids in porous media. One of the main objective is to derive an generalized Darcy equation for yield stress fluids relating mean flow rate, the pressure gradient and a critical pressure gradient. Experimentally [2, 9, 19] they establish various laws which have quite essentially the same structure,
where Q is the total flow rate and n is the rheological exponent of the Herschel-Bulkley model. In particular, n = 1 defines the Bingham fluid. However, numerical simulations and pore-network models [15, 17, 20] demonstrate the situation may be much more complex that Equation (1) may be expected when the flow occurs in non-trivial geometries. Three different flow regimes,
with α successively equal to 1, 2, and 1 for Bingham fluids when ΔP is increased beyond ΔPc, may be identified. These three regimes are a consequence of the presence of heterogeneities in the pore size distribution and network topology.
In the present paper, we will study the different flow regimes in a one-dimensional (1D) rough channel. Besides offering interesting insights into the complexitites of the flow of yield stress fluids even in this seemingly simple case, it is of direct relevance in connection with flow both in two and three-dimensional porous media flow. Close to the yield threshold, there will be only one single one-dimensional channel where there is flow irrespecive of the dimensionality of the porous medium (see e.g., [11, 15, 17, 21]).
The present work relates to earlier work considering channel flow with periodically varying aperture [22–24]. In those articles, the authors investigated the validity of the lubrication approximation  which can become inconsistent with the presence of a unyielded plug along the center of the channel.
The main purpose of this paper is to investigate how Bingham fluids behave in one-dimensional channels with self-affine aperture variations. From earlier studies of how the permeability is affected by self-affine correlations in the aperture field , non-trivial scaling is to be expected of the parameters that characterize the flow also in this case.
A channel with a self-affine aperture field h(x) has an average aperture
where L is the distance along the channel over which the average is taken and H is the Hurst exponent . We will later on this this paper present a more precise definition. Fractures are typically self affine (see e.g., [27–30]).
In the case of Newtonian flow in a channel with self-affine aperture correlations, it was enough to analyze the scaling properties of the permeability to have a full overview of what effect such correlations have on the flow. In the present case, it is necessary to expand the analytical solutions in power series. Each term in this series shows different scaling behavior.
The outline of the paper is as follow. In section 2, we consider flow in a one-dimensional rough channel in the lubrication limit. We solve the flow equations analytically and find flow rate as function of pressure difference across the length of the channel. After identifying a length scale associated to the flow rate, we use this length to expand in power series for small and large values of it, our analytical solution. In section 3, we discuss briefly our solution for a channel aperture that follows a cosine variation. In section 4, which contains the main results of this paper, we consider the roughness to be self affine and derive scaling relations between aperture, length of the channel and the length scale associated with the flow rate. In comparison to the Newtonian case, the scaling properties are much more subtle in this case due to the presence of the length scale, a*, characterizing the flow. We use the series expansions of section 2 to determine how each term scales. Finally, we propose in section 5 a truncated version of the double series expansion which reproduces the analytical solution over the entire range of flow rates. This is very useful as it is necessary to use such series to characterize the scaling properties of the flow equations when the aperture field is self affine. Section 6 contains our conclusions.
2. Bingham Flow in a Rough Channel: General Results
where xy is the shear rate, σxy is the shear stress, σc is the shear stress yield threshold, D is the viscosity and sgn is the sign function.
By integrating the flow equations across the channel with the constitutive law given in Equation (4), we find for for position x along the channel where the aperture is h = h(x) that the total flow rate Q per unit of width is given by
where p′ = p′(x) is the pressure gradient at x. We have also defined a local pressure gradient threshold p′c = p′c(x) given by
We have here implicitly made the assumption that h(x) varies slowly enough so that the lubrication limit is in effect .
We may non-dimensionalize the flow Equation (5) by introducing a dimensionless pressure
and the inverse of the local Bingham number:
Which quantify the balance between the viscous and yield stresses contributions at the scale h. Equation (5) then reads
We solve this equation with respect to p′ and find
By integrating this equation over the length of the channel, we express the total flow rate Q as a function of the pressure drop over the total channel length L,
Hence, using Equations (7), (8), and (10) we find
We note that |ΔP| → Pc as |Q| → 0 in this equation, where
We have here used Equation (6) and defined
Pc is the minimum pressure drop needed so that |Q| > 0 when |ΔP| > Pc. It is worthwhile to derive this threshold pressure from direct physical arguments since the result is not obvious.
Let us pick two positions along the channel, one at x1 and the other at x2. Here the threshold pressures are respectively p′c1 and p′c2. Since this is a channel, the total flux Q is the same everywhere along it. With reference to Figure 1, we now assume that the total flux is Qa. Correspondingly the pressure gradient at position xa is p′1a and at position x2 it is p′2a. If we lower the total flux to Qb we respectively find corresponding lower pressure gradients p′1b and p′2b at x1 and x2 respectively. We see from the figure that as the total flux Q → 0, p′1 reach the threshold value p′c1Âăexactly as p′2 reaches p′c2. Hence, as the total flux reaches zero all along the channel, the local pressure gradients reach their threshold values simultaneously. Equation (13) then follows since the total pressure as Q → 0 then must be the integral over the local threshold pressure gradients.
Figure 1. Labels “1” and “2” refer to two positions along a one-dimensional channel. The pressure gradient thresholds are p′c1 and p′c2 respectively. The total flow rate through the channel is the same at all positions. Hence, when the total flow rate is Qa, there will be a corresponding pressure gradient p′1a at position 1 and a pressure gradient p′2b at position 2. When the total flow rate is lowered to Qb, the pressure gradients become p′1b and p′2b at positions 1 and 2 respectively. From the figure, we see that as Q → 0, the pressure gradients p′1 and p′2 reach their respective thresholds, p′c1 and p′c2, simultaneously.
We now return to Equation (10). We may expand it for small |Q|. We find
where a0 = 1, a1 = , a2 = 16/3, and a3 = 40/9.
If we keep only two terms in this series, reintroduce the original variables and integrate along the x axis, we find
and Pc is given by Equation (13). Hence, close to the threshold pressure, the flux is quadratic in the excess pressure.
Let us now introduce a length scale a* defined as
We have used here that Q is a constant along the length of the tube, and therefore Q varies with x. a* characterizes the flow everywhere along it and is independent of x. Hence, |Q| ≪ 1 is equivalent to a* ≪ h(x) and vice versa |Q| ≫ 1 is equivalent to a* ≫ h(x) for all x.
The series expansion (15) may then be transformed into the small-a* expansion
where we have defined
where we have split the channel aperture into two components,
All the geometrical information on the channel is now found in the integrals Ik(a, L). The information on the flow is encoded in the length scale a*.
In the other limit, Q ≫ 1, Equation (10) may be expanded as
where b2 = −1/288, b3 = 1/1152 and so on. We introduce the length scale a* (18) in this expression and find
which is then a large-a* expansion.
With only the dominating term, proportional to (a*)2, included in this equation, we recover the Darcy equation with viscosity ν = D. Hence, the permeability is controlled by the integral I3(L), as has already been discussed in Balhoff et al. .
It is interesting to note that the term corresponding to Darcy flow is present both in the large a* and the small a* expansions. It needs to be present in the large-a* expansion even though it is proportional to (a*)2 since the Bingham fluid approaches the Newtonian limit with increasing flow rate.
Keeping the two leading terms and reverting to the original variables gives the equation
where Pc is the threshold pressure. Hence, the fluid behaves as a linear fluid, but with an effective pressure threshold that must be overcome for flow to occur. The effective threshold is larger than the real flow threshold Pc as one would expect. Equation (26) is the approximation used e.g., by Roux and Herrmann  in their study of Bingham fluids in disordered networks. Since the effective pressure threshold, Equation (27) is proportional to the real threshold, it is still given by Equation (14). This justifies the approach of Roux and Herrmann .
3. Periodic Aperture Variations
Both expansion (19) and (25) are applicable for any aperture distribution, as long as the lubrication approximation remains valid. As example, following Frigaard et al. [23, 24] who considered an aperture varying according to a cosine function. We may then calculate Ik(a, L) analytically. We consider k = 1, 2, and 3.
Considering the aperture
where b is the amplitude of the aperture variation and N is the number of wavelengths in all. We then have
Hence, the threshold pressure is then
and the two limit expansions (19) and (25) to leading order are
4. Self-Affine Apertures
A large number of studies report that the roughness of fracture surfaces show self-affine correlation [27–30]. We consider here the flow in a channel with aperture h(x) with a self-affine properties. This means that the aperture is characterized by a two point function p2(h2 − h1, x2 − x1) giving the probability density that if h(x1) = h1 then h(x2) = h2. Self affinity is the scaling property
Normally, we have that 0 < H < 1.
In an earlier study , we considered the flow of a newtonian fluid in a channel where the aperture field obeyes the scaling relation (33). After identifying a length scale a characterizing the aperture opening (to be defined below), the permeability K relating flow rate Q to pressure difference along the channel, ΔP, scales as K ~ aκ. For small a, κ = 3, for intermediate a, κ = 3 − 1/H and for large a, κ = 3. The question we pose here is: how does the the self affinity of the aperture field influence the equation relating Q and ΔP for a Bingham fluid.
The channel aperture is split into two components, h(x) = a + η(x), where the minimum aperture is a = minx h(x) and η(x) is self affine characterized by a Hurst exponent H. Since we are in the lubrication limit, we may reshuffle the aperture field η(x) → η[ξ] = η(x(ξ)) where we have used an ordering transformation such that η(ξ1) ≤ η(ξ2) if ξ1 < ξ2. The averaged ordered sequence obeys the scaling law 
where we have set a prefactor equal to unity.
The fundamental integrals (20) that control the transport equations may then be written
where 2F1 is the Gaussian hypergeometric function.
First we assume k > 0 in (35). For large a ≫ LH, the integral is proportional to La−k.
For intermediate to small values of a ≪ LH, we may rewrite the integral
Depending on whether kH is smaller than or larger than one, the integral will be dominated by the upper or lower integration limit. Hence, if kH < 1, the integral behaves as
This leads to
On the other hand, if kH > 1, the integral will be dominated by the lower limit, and it does not depend on L nor a. We may summarize this,
When a becomes so small that the region around the bottle neck (where the aperture is a) dominates, integral (36) must be discretized ,
where Δ is the discretization. For a << (Δ/L)H, this integral becomes dominated by the j = 0 term and we have
We may summarize this discussion through the equation
The behavior of Ik(a, L) for negative k is simpler. We only have two cases to consider: a ≫ LH and a ≪ LH. We find
We test in Figure 2 the scaling in Equation (42).
Figure 2. Top: Numerical computation of integrals Ik(a, L) defined in Equation (20) for two different Hurst exponents (plain line: H = 0.8, dashed line: H = 0.2). The size is L = 219 and the integral has been sample averaged over twenty realizations. The regime of interest, where non-trivial scaling will occur, is in the region of intermediate values of a. As the figure runs over 11 orders of magnitude in a, this region may seem narrow in the figures. Bottom: The integrals have been rescaled by akH−1 Ik(a, L) in order to bring out the scaling regions more clearly. The region where we see scaling according to Equation (42) now appears horizontal. We have chosen H = 0.8.
In the three upper panels, we show Ik(a, L) with k = 1, 2, and 3 for H = 0.2 and H = 0.8. For H = 0.2, kH is always less than 1, the second regime is thus independent of a and the top curves of Figure 2 for H = 0.2 (dotted lines) show a flat middle section. For H = 0.8, we have kH < 1 only for k = 1 whereas for k = 2 and 3, we have kH > 1, the two curves for H = 0.8 and k = 2 and 3 should not, thus, show a flat middle section; this is actually what is observed in Figure 2. In the lower three panels we show the rescaled integrals akH − 1 Ik(a, L) for k = 1, 2, and 3. Equation (42) then predicts that the middle sections of the two curves corresponding to H = 0.8 k = 2 and 3 (values for which kH > 1) should be flat whereas those for H = 0.2 and (H = 0.8, k = 1) should not (in those cases we have kH < 1).
We now return to the series expansions of the relation between pressure drop ΔP and flow rate Q. We have introduced the length scale a* in Equation (18) and the expansion for small flow rates—hence, small a*—is given in Equation (19). Combining this expression with the assumption that ΔH ≪ a ≪ LH, we find
where [·] is the floor function. Hence, for 1/2 < H < 1, only the k = 0 term scales with L. For 1/3 < H < 1/2, the k = 0 and k = 1 terms scale with L.
The k = 0 term in (44) is the threshold pressure. We therefore have that
The first regime corresponds to the situation where the yield stress is controlled by the minimal aperture region. In that case, the pressure gradient is dominated by the pressure drop over this particular region of the length Δ. The second threshold corresponds to the limit where yield stress is controlled the full length of the rough channel; in that case the typical aperture is related the self-affine nature of aperture and scales with the length of the channel like LH. The third asymptotic regime is a channel of length L with parallel walls separate by a distance a. Let us now consider the large |Q|—the large a*—expansion. By combining Equation (25) with the scaling of Ik(a, L), we find
as long as H > 1/3.
The scaling seen in this section is subtle. The difficulty lies in the presence of the length scale a*. This length scale is not present when the flow is Newtonian. In that case, the expression for the pressure drop vs. flow rate is 1
where μ the viscosity. Hence, there are three scaling regimes in this case .
5. A Truncated Series Expansion
We present in this section a truncation of the series expansions of the flow equation, which nevertheless retains very close to the exact expression over the entire range of aperture opening a. This series is considerably simpler than the exact solution.
We combine the two expansions (19) and (25) by splitting up the integrals Ik(a, L), defined in Equation (20),
We combine the expansions (19) and (25) and get
We truncate the series at second order, finding
where we have used I1(< a*, L) + I1(> a*, L) = I1.
In Figure 3, we plot this approximation and compare it to the exact solution.
Figure 3. Flow rate Q as function of the applied pressure drop ΔP − Pc for different Hurst coefficient H and different openings a. The other parameters are D = 1/12, σc = 1, L = 219. Black plain line represent numerical integration of the exact expression of Equation (12). Dashed lines represents the asymptotic regimes for Q ≪ 1 defined by Equation (44) (thick line: second order for thin line for first order) and for Q ≫ 1 defined by Equation (46) (first order). Open circles represent the combined asymptotic approximation Equation (50). The horizontal line, represents the flow rate at which a* = a, which represents the first decrease (resp. increase) of Ik(< a*, L) (resp. Ik(> a*, L). This flow rate corresponds also to the maximal error between the exact and approximated solution. This errors is respectively (from top-left to bottom-right): 10, 2, 12, and 5%.
The agreement with the analytical summation is very good for any opening and any Hurst coefficient. The reason for this good match can be understood as follows. At very low flow rate, if a* is smaller than any openings, we have then Ik(< a*, L) = 0 and Ik(> a*, L) = Ik(a, L). The combined expression is exactly equal to the expansion of Equation (44). As Q, one start to reach the point where a* = a (denoted as an horizontal line in Figure 3) which corresponds to the situation where the smallest element of the sum Ik(> a*, L) is removed but added in Ik(< a*). From there, all the integrals Ik(> a*, L) (respectively Ik(< a*, L)) are decreasing (respectively increasing) with Q until a* reaches the maximum opening. For this case, one then have Ik(> a*, L) = 0 and Ik(< a*, L) = Ik(a, L) and we recover the limiting case of Equation (46). The combined approximation allows us thus to match the two asymptotic limits of Equation (44) and (46).
This good agreement between the approximation and the exact solution let us to deduce sufficient condition of validity for the two expansions. Namely:
The validity of the regime reflects thus the heterogeneities of the aperture, particularly the two extremal values.
In this work, we have investigated the flow rate in a 1D rough channel for any aperture distribution. We have, however, paid particular attention to the case of aperture fields that are self affine, which is e.g., the case when the aperture field is given by the channel being a fracture.
Our first result is given in Equation (12). This gives the total flow rate Q as a function of pressure difference ΔP across the length of the channel. Whatever the geometry of the channel, close to the pressure threshold Pc, the flow rate scales like (ΔP − Pc)2 with a prefactor function of the fluid and aperture properties. We have introduced a characteristic length scale a* = where the square root scaling is a direct consequence of the square variation of the flow rate with the pressure. The length scales a* allows an expansion in series of the two limit cases, the low (a* < a) and high (a* > a) flow rates regimes where a is the minimal aperture. At large flow rate (large a*), the flow rate scales linearly as (ΔP − Pc), with an apparent critical pressure Pc larger than the actual one Pc.
With this length scale, we combine the series expansions and truncate them in a way that leads to high accuracy. The sum of the two contributions is found to provide a good approximation for the exact solution given by Equation (12) in the case of self-affine channels.
One of the main consequences of this results is the fact that, in principle, one should expect four different flowing regimes in porous media. Presently, only three flowing regimes have been observed numerically (see ). The first one, in the limit ΔP ~ Pc, is linear ( i.e., Q ~ (ΔP − Pc)) and corresponds the flow channelization along one channel. Here we have demonstrate that, in addition to this regime, one should observed another quadratic regime where Q ~ (ΔP − Pc)2. As a consequence, by increasing the pressure, one should observes regimes in the following sequence: quadratic (one channel), linear (one channel), quadratic (with an increase of channels) and linear (when all channels are open).
It is important to recall that this work is only valid under the lubrication limit. Lubrication approximations have been investigated extensively for Newtonian fluids in confined geometries (e.g., [32, 33]). For visco-plastic fluids, the limit of lubrication has also been investigated by various authors in different geometries (see  for film flows, [22–24, 35] for confined geometries and [36, 37] for experiments). In this case, the lubrication hypothesis becomes very restrictive essentially due to rigidity of the plug regions. As discussed above (and in details by Frigaard and Ryan ), because of the aperture variations and the mass conservation, lubrication approximation leads to variations of the plug velocity. However, such variations are forbidden in the case of ideal visco-plastic fluids which prevent any extensional deformations. Consequently, one expects deviations from the lubrication hypothesis as the amplitude of the variations becomes important.
As for flow in channels with self-affine aperture fields, we do not find a simple scaling law for the equivalent of the permeability as was the case for newtonian fluids . Rather, we find the each term in the series expansions of the flow equation scales with its own scaling exponent in the aperture opening a. By combining these results with the truncated series given in section 5, it is tractable how the flow equations evolve with varying a in this case.
Future work will be dedicated to the investigation of the validity conditions of our analytical expressions. We have not considered two and three-dimensional channels. As for the case of Newtonian flow, this is much more complex than the present case. However, as was remarked in the Introduction, at small pressure gradients, the flow will be confined to single channels, i.e., the precise situation considered in the present paper.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors would like to thank J. P. Hulin, S. Roux, and D. Salin for interesting for useful discussions, the “Agence National de la Recherche” for financial support of the project LaboCothep ANR-12-MONU-0011.
2. Chevalier T, Chevalier C, Clain X, Dupla JC, Canou J, Rodts S, et al. Darcy's law for yield stress fluid flowing through a porous medium. J Non-Newtonian Fluid Mech. (2013) 195:57–66. doi: 10.1016/j.jnnfm.2012.12.005
15. Talon L, Bauer D. On the determination of a generalized Darcy equation for yield-stress fluid in porous media using a Lattice-Boltzmann TRT scheme. Eur Phys J E (2013) 36:1–10. doi: 10.1140/epje/i2013-13139-3
16. Balhoff M, Sanchez-Rivera D, Kwok A, Mehmani Y, Prodanović, M. Numerical algorithms for network modeling of yield stress and other non-Newtonian fluids in porous media. Trans. Porous Media (2012) 93:363–79.
22. Putz A, Frigaard IA, Martinez DM. On the lubrication paradox and the use of regularisation methods for lubrication flows. J Non-Newtonian Fluid Mech. (2009) 163:62–77. doi: 10.1016/j.jnnfm.2009.06.006
24. Roustaei A, Frigaard IA. The occurrence of fouling layers in the flow of a yield stress fluid along a wavy-walled channel. J Non-Newtonian Fluid Mech. (2013) 198:109–24. doi: 10.1016/j.jnnfm.2013.03.005
32. Reynolds O. On the theory of lubrication and its application to Mr. beauchamp tower's experiments, including an experimental determination of the viscosity of olive oil. Philos Trans R Soc Lond B Biol Sci. (1886) 177:157–234.
36. de Souza Mendes PR, Naccache MF, Varges PR, Marchesini FH. Flow of viscoplastic liquids through axisymmetric expansions contractions. J Non-Newtonian Fluid Mech. (2007) 142:207–17. doi: 10.1016/j.jnnfm.2006.09.007
Keywords: non-newtonian flow, yield-stress fluid, rheology, flow in fracture, lubrication
Citation: Talon L, Auradou H and Hansen A (2014) Effective rheology of Bingham fluids in a rough channel. Front. Physics 2:24. doi: 10.3389/fphy.2014.00024
Received: 16 December 2013; Accepted: 02 April 2014;
Published online: 23 April 2014.
Edited by:Colm O'Dwyer, University College Cork, Ireland
Reviewed by:Markus Wilhelm Abel, Ambrosys GmbH, Germany
Allbens Picardi Faria Atman, Centro Federal de Educação Tecnológica de Minas Gerais, Brazil
Copyright © 2014 Talon, Auradou and Hansen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Laurent Talon, CNRS, FAST, Université de Paris-Sud, Bât. 502, Campus University, Orsay, F–91405, France e-mail: firstname.lastname@example.org