The Effect of Microstructure on Models for the Flow of a Bingham Fluid in Porous Media: One-Dimensional Flows

The Buckingham–Reiner models for the one-dimensional flow of a Bingham fluid along a uniform pipe or channel are well-known, but are modified here to cover much more general one-dimensional configurations. These include selections of channels with different widths, and five different probability density functions describing distributions of channel widths. It is found that the manner in which breakthrough occurs at the threshold pressure gradient depends very strongly on the type of distribution of pores and that a pseudo-threshold pressure gradient, which might be inferred from measurements of flow at relatively high pressure gradients, may be more than twice the magnitude of the true threshold gradient.


Introduction
Many fluids which are used in everyday life and in industry exhibit a yield threshold, by which is meant that the fluid either remains stationary, or else moves as a solid without shear, whenever the shear stress is lower than a threshold value. Fluids which exhibit such a threshold include Casson fluids, Herschel-Bulkley fluids and Bingham fluids. Of these, the Bingham fluid exhibits a linear stress-strain relationship post-threshold. While the literature which is associated with such non-Newtonian fluids is extensive and mature, the research associated with the flow of such fluids through porous media is in the most part relatively small and recent. This short paper, then, is devoted to the task of determining how Darcy's law, which governs the flow of a Newtonian fluid in a porous medium, is altered when the porous medium is saturated by a Bingham fluid. Many authors have considered the porous medium to consist either of bundles of tubes or networks of tubes of circular cross section and of various diameters. Flow within such tubes was one of the main considerations in the lengthy work of Bingham (1916), which describes the experimental investigations by Simonis (1905) on the rate of flow of suspensions of fine clay in NaOH solutions through tubes of different radii. Bingham concluded that these suspensions have a yield threshold (i.e. an applied pressure gradient above which the fluid flows) which he called the friction constant, and he stated that the response to the pressure is then linear post-yield. Such a simple piecewise-linear model of the flow dependence on the applied pressure gradient was then assumed to apply for flows in porous media in work by Pascal (1979Pascal ( , 1981Pascal ( , 1983. Bernardiner and Protopapas (1994), in a review of pre-1994 Russian literature, also assume this model; they cite experimental work involving viscoplastic fluids, and Newtonian flows through argillaceous soils, the latter being very similar to Bingham's suspensions.
Despite Bingham's stated conclusion that the post-yield flow rate depends linearly on the applied pressure gradient, his Fig. 12 (which was taken from Simonis 1905 and 'adapted' by Bingham) shows a well-defined curved section immediately post-yield; this, he concluded, was due to seepage, suggesting implicitly that its presence was an experimental error. In slightly later, and apparently independent papers, Buckingham (1921) and Reiner (1926) derived the now well-known and more recently named Buckingham-Reiner law for flow in cylindrical tubes. This formula, seen later, matches well the shape of Simonis's curves at least qualitatively, therefore, it might be concluded that Bingham's supposed reason for the experimental error in the data of Simonis no longer needs to be invoked.
Others have considered different types of channel flow and different aspects of those flows. For example, Daprà and Scarpi (2004) considered unsteady flow in a plane channel, which might be termed a Poiseuille-Bingham flow, and developed a series solution for the start-up flow from rest using the Pascal model. This was extended to pipe flow in Daprà and Scarpi (2005). Muraleva and Muraleva (2011) considered steady flow in an undulating channel using finite difference methods to determine the unyielded regions. Talon et al. (2014) considered more general boundary imperfections for the channel within the lubrication approximation, and they applied their method to periodic and self-affine undulations. The computation of pipe flow in more complicated geometries, such as square cross sections, is hindered by the need to find the yield surfaces numerically. Considerable progress has been made on this aspect, and we refer the reader to the work of Mosolov and Miasnikov (1966), Saramito and Roquet (2001), Huilgol (2006) and Muraleva and Muraleva (2009). It is worth mentioning that Saramito and Roquet (2001) show that the Buckingham-Reiner law for a circular tube may be used with minimal error for a square tube after rescaling. On the other hand, pore networks have also been used to simulate the flow of Bingham fluids in porous media. The paper by Balhoff et al. (2012) reports on some different numerical methods for computing the flow in networks. Flows within a diamond-shaped network were considered by Chevalier and Talon (2015) where the different paths were assigned random thresholds. Attention was then focussed on the numerical simulation on the successive initiation of flow within the different paths from one side of the porous medium to the other. Another analysis of one-dimensional flow within a two-dimensional network was considered by Chen et al. (2005), but this was in the context of imbibition and the displacement of a Bingham fluid by a Newtonian fluid. We would also wish to draw the readers' attention to the papers by Chevalier et al. (2013Chevalier et al. ( , 2014, which report on the experimental determination of a suitable Darcy's law for Herschel-Bulkley fluids, and the paper by Bleyer and Coussot, which is concerned with the numerical modelling of a Bingham fluid through an array of circular cylinders. Entrance effects have also been considered by Philippou et al. (2016).
The present paper has multiple motivations. One of these is the experimental work of Chase and Dachavijit (2005), which uses a packed column containing spherical glass beads, where it is found that a scaled Buckingham-Reiner equation may be used to model at least this type of particulate porous medium. A second is the analysis of Oukhlev et al. (2014) which shows how a Bingham fluid may be used to determine the pore size distribution within a medium. Their method was tested successfully on a medium consisting of pores which satisfy a Gaussian distribution. An analogous paper by Castro et al. (2014) seeks to determine the pore size distribution by comparing the increased flow rate caused by successive increments in the applied pressure gradient. A third might seem negative, but Coussot (2014) states that network models (a tube bundle being a subset of these) do not provide generic analytical expressions for the velocity/pressure-drop relationship. The present study will show analytically that this relationship does indeed depend quite strongly on the type of distribution of channels/tubes which are present, and therefore Coussot's view will be confirmed.
Given that this is the centenary of Bingham's original paper (at the time of submission of the present manuscript, that is), it seems apposite, then, to revisit the simplicity of the onedimensional channel/tube model of a porous medium in order to show analytically how the Buckingham-Reiner model depends very strongly on the distribution of the channels/tubes. In the main body of the paper, we consider various distributions of channels, both discrete and continuous, while the "Appendix" section lists the equivalent formulae for circular tubes. We find that the manner in which the mean velocity varies with the applied pressure gradient immediately post-threshold takes various power-law forms depending on the distribution of channels/tubes. It is also found that the quantity which we shall term the pseudo-threshold (i.e. the threshold gradient which is given when extrapolating back from flows at high pressure gradients) also varies greatly with the type of distribution of channels/tubes.

The plane-Poiseuille flow of a Bingham fluid
A Bingham fluid is characterised by having a yield stress, but it is otherwise Newtonian, and therefore it flows only when a sufficiently large stress has been imposed. This stress may arise from pressure gradients, moving boundaries or buoyancy forces, for example. The simplest case in which this happens is for Poiseuille flow where we shall take u(y) to be the velocity in the x-direction. For one-dimensional isothermal flow in the x-direction, we have where τ is the shear stress. If τ 0 is the yield stress, then the rate of strain is given by In what follows, we shall work with G = −d p/dx, which is the pressure drop per unit length, and therefore a positive value of G causes the velocity to be positive whenever it is nonzero. If a uniform channel of width, h, has bounding planes at y = ±h/2, then the velocity profile in the channel takes the form where the value, , is given by and where this solution is valid only when < 1. The value of is the proportion of the channel which is unyielded, and therefore values which are above unity mean that the whole channel is unyielded and hence u = 0. It proves to be much more convenient to work with because this may be interpreted as a nondimensional pressure gradient, one which is based on the yield stress of the fluid. Thus, the fluid begins to flow once σ rises above unity in magnitude, and so σ = 1 is the nondimensional threshold pressure gradient. The velocity profile given by Eq. (3) is parabolic in two outer regions and is a nonzero constant in the middle region where the unyielded fluid forms a plug.
In the context of porous media, it is necessary first to find the mean flow. Therefore, we shall assume in the first instance that the porous medium consists of a periodic array of these channels where the period is H . The porosity of such a medium is then φ = h/H . Within one channel the total volumetric flux, Q, per unit distance in the third coordinate direction (z) is given by this defines the function f (σ ) which is the plane-Poiseuille flow equivalent of the widely known Buckingham-Reiner formula for the Hagen-Poiseuille flow of a Bingham fluid; see the "Appendix" section. In the remainder of this paper, we will omit the use of the modulus signs in formulae such as Eq. (6), for simplicity of presentation, and will restrict attention to positive values of G and hence of σ . Given the formula in Eq. (6), the equivalent Darcy velocity for a Bingham fluid, u B , is the mean velocity taken over one period of the pattern of channels. Hence, where are the permeability and porosity, respectively, of the present porous medium. In what follows, the definitions of K and φ will depend strongly on the choice of channels which are used to define one period of the porous medium. The general guiding principle will be that the σ 1 limit, which corresponds to Newtonian flow, should always take the form, u N = G K /μ. This means that we have the following comparison between the response of Bingham and Newtonian fluids for the present configuration: The variation of f (σ ) is shown in Fig. 1a. While its nature is such that f (σ ) → 1 as σ → ∞, its behaviour near to the threshold takes a quadratic form at leading order, Therefore, the induced velocity first begins to rise from zero in a quadratic manner once the threshold pressure gradient is exceeded; this is precisely the qualitative behaviour found in Fig. 12 of Bingham (1916). At the opposite extreme of large values of σ , a Newtonian behaviour is expected. However, the induced flow is equal to σ f (σ ) in nondimensional terms (see Eq. (6)) and this is shown in Fig 1b. The two-term large-σ approximation to this flow is given by, From an experimental point of view, if one were to try to extrapolate back from fairly large values of σ to try to find the threshold gradient, then Eq. (11) suggests that the value σ = 3/2 might be found because σ f (σ ) tends quite quickly to the linear behaviour shown in Eq. (11) as σ increases above 1; see the dotted line in Fig 1b. We may therefore define a pseudo-threshold value of for a porous medium consisting of identical channels; this value may also be found in Talon et al. (2014) where it is termed the effective pressure threshold or apparent critical pressure. The equivalent analysis for Hagen-Poiseuille flow is summarised in the "Appendix" section where the analogous function g(σ ) is defined and which applies to a circular pore of radius, h. This curve is shown in Fig. 1 as a dashed line; it displays very similar characteristics to that of f (σ ) with an initial quadratic rise as σ increases above 1 and where it tends towards 1 as σ → ∞.
In the rest of the paper, we will use F(σ ) to denote u B /u N , the ratio between the flow rate of a Bingham fluid and that of a Newtonian fluid, and therefore the nondimensional Darcy flow is σ F(σ ). The equivalent notations for Hagen-Poiseuille flow in circular channels, which is summarised in the "Appendix" section, are G(σ ) and σ G(σ ).

The effect of having two different channel widths
A porous medium will generally be formed of channels with different widths, unless they have been designed otherwise, and therefore we begin a generalisation of the above analysis of the type of Purcell (1949) who considered a porous medium consisting of cylindrical pores with many different diameters. We will now assume that one period of channels will consist of one channel of width, h, and a second narrower channel of width γ h where γ < 1. Clearly, if σ = 1 represents the threshold pressure gradient for the wider channel, then σ = 1/γ is the threshold for the narrower channel. Therefore, the total fluid flux through one of each channel is, The corresponding Darcy velocity is now given by, or where the permeability and porosity are given by The pseudo-threshold value of σ is Samples of the flow relative to that of a Newtonian fluid are shown in Fig. 2. Attention is focussed on γ = 0.25, 0.5 and 0.75, and the single-channel case, which may be regarded as being perfectly equivalent to either γ = 0 or γ = 1, is included as a reference case. When σ takes relatively small values, F(σ ) appears to decrease as γ increases from zerothis may be expected since flow arises in only the wider channel and hence Eq. (15) yields F(σ ) = f (σ )/(1 + γ 2 )-but eventually flow begins in the second channel, we return to the single-channel case as γ → 1, and clearly the variation in F(σ ) with γ cannot be monotonic. Indeed, the F(σ ) curves for γ = 0.5 and γ = 0.75 close to σ = 3.2. Figure 3 allows us to see more clearly why there is this lack of monotonicity in the variation of F(σ ) with γ . Each curve corresponds to a chosen value of σ between 1.5 and 100, and each has two symbols. The first (o) shows the value of γ above which the narrower channel admits flow. The second (•) depicts the smallest value of F(σ ) as γ varies. It is to be expected that, when the permeability is correctly accounted for, the value of F(σ )n must be identical in the two cases, γ = 0 and γ = 1, because they describe effectively the same situation, namely that all channels have width, h. We have already seen that the (1 + γ 3 ) scaling factor causes F(σ ) to decrease as γ increases when flow arises only in the wider channels. Once flow begins in the narrower channels, this decrease is attenuated, and given that flow in the narrower channels becomes easier as γ increases, F(σ ) must then achieve a minimum before γ eventually rises up to unity. Figure 3 also suggests that the value of γ at which F(σ ) has its minimum tends towards 1 as σ → 1 + . This may be analysed by setting σ = 1 + δ σ and γ = 1 − δ γ into Eq. (15) and expanding for small values of both δ σ and δ γ , both values being assumed to be positive. At leading order, we find that which shows that, for a fixed value of δ σ , the minimum value of F(σ ) occurs when δ γ = δ σ as σ → 1 + . For Hagen-Poiseuille flow, the equivalent equation is which was obtained from Eq. (53), and therefore exactly the same conclusion applies.

The effect of having many different channel widths
The above two-channel configuration may be generalised further by considering N channels per period. If the widths of the channels are taken to be γ i h for i = 1, . . . , N , and where the channels have been ordered with respect to their widths, then we have This definition allows one to define simple ratios of channel widths. For example, if one were to take five channels per period where three have width, h, and the remaining two have width, h/2, then we may set γ 1 = γ 2 = γ 3 = 1 and γ 4 = γ 5 = 1 2 with N = 5. The Darcy velocity for the general case is now given by while the permeability and porosity are given by, The pseudo-threshold is Some sample cases of multiple channels are given in Figs. 4 and 5. Figure 4 is concerned with the effect of having different numbers of channels of width, 1 2 h, in addition to one channel of width, h. We see that the initial curves in the range 1 < σ < 2 differ from one another and are dependent on the number of smaller channels, despite flow taking place only in the widest channel. In practice, the overall mass flux will be identical, but we are plotting the flow which is obtained relative to that of a Newtonian fluid for which there will be flow in the narrower channels. Mathematically, the reduction in F(σ ) is due to the presence of the denominator in Eq. (21). The threshold pressure gradient in the narrower channels is σ = 2, and this only becomes very clearly apparent in the behaviour of the curves when there are a large number of narrower channels; see, for example, when n = 16 and n = 25. Figure 5 shows a selection of cases which approximate a uniform distribution of channels, a case which will be considered in the next section. The general configuration consists of n channels where the values of γ take the form, γ i = i/n, for i = 1, . . . , n. We show curves which correspond to n = 5, 10, 20 and 40, and the uniform distribution response given by Eq. (30). The single-channel case, which could be said to be equivalent to n = 0, is given in Fig. 5a for reference.
In general, the larger the value of n, the smaller is the value of F(σ ), but the curve representing the variation of F(σ ) tends monotonically towards that of the uniform distribution, Eq. (30); see Fig. 5a. In particular, the initial (σ − 1) 2 dependence of F(σ ) which is given in Eq. (10) weakens as n increases as only one out of many channels admits flow followed by another and then others as the pressure gradient increases. The limit as n → ∞, namely the uniform distribution, then represents an initial cubic dependence on (σ − 1), as will be shown below. Figure 5b shows that, even for a discrete set of 20 channels with a uniformly distributed set of widths, a continuous uniform distribution of channels may be modelled closely.

Random distributions of channels
More likely again is that the porous medium should consist of a distribution of channel widths, and while one may investigate many different types, we shall confine ourselves solely to five, namely a uniform distribution, two different types of linear distributions, a quadratic distribution and a seminormal distribution. We may use the function, Φ(γ ), as the probability density function for the distribution. In general, expressions for the Darcy velocity relative to that of a Newtonian fluid, the permeability and the porosity take the following forms, where

Uniform distribution
For a uniform distribution of pores in the range, a ≤ γ ≤ 1, the probability density function is, therefore, Guided by the above analyses, the Darcy velocity for such a distribution of pores is For this case, the permeability and porosity are given by and the pseudo-threshold is .
Thus, there are three regimes depending on the value of σ . The first region is stagnant, while the second (intermediate) regime has flow arising in a proportion of the channels (i.e. those with widths lying between h/σ and h), and then the final regime corresponds to flow in all channels. Figure 6 shows how F(σ ) varies with σ for a selection of values of a. In all the cases where a = 1, it is important to note that flow starts very weakly as σ increases above unity. The formula in Eq. (27) for the intermediate regime shows that the initial increase in the velocity is now proportional to (σ − 1) 3 , as opposed to the single-channel case where the increase is quadratic. We also note that, when a = 0, i.e. when channel widths do not have a lower limit, then the final regime does not exist. In this case, we have the relatively simple relation In this limit, the pseudo-threshold becomes σ pt = 2. This case is difficult to see in Fig. 6 because the curve is almost indistinguishable from that of a = 0.25. On the other hand, when we formally let a → 1, which corresponds to almost no variation in the channel widths, the intermediate regime no longer exists, and the formula corresponding to the final regime in Eq. (27) reduces precisely to that given in Eq. (9) for the single channel. Both of these limits are seen in Fig. 6. Similarly, the location of the pseudo-threshold (Eq. (29)) reduces to that for a identical channels; see Eq. (12).

Linear distributions
We now consider the following two linear distributions of channel widths for which the probability density functions are, The first corresponds to situations where very narrow channels are highly unlikely, while the second makes the narrow channels very likely. For case (a), our analysis gives where the permeability and porosity are given by, The pseudo-threshold is For case (b), we have where K = 1 10 In each of the two linear cases, we recover a unit scaled velocity as σ → ∞, but the onset of flow near σ = 1 is marked by a different power of (σ − 1). Thus, for case (b), flow In Frame a, the curves correspond to the following distributions: single channel (C1); uniform distribution with a = 0 (U ); the two linear distributions (L1 and L2), the quadratic distribution (Q) and the seminormal distribution (SN ). In the three remaining frames, all the curves are in the same order. Also shown as dotted lines are the extrapolations leading to the locations of the pseudo-thresholds begins very weakly indeed because of the additional lack of availability of channels of width h within which to flow as compared to the uniform distribution case. This also results in an increased value for the pseudo-threshold: Both of these conclusions may be seen clearly in Fig. 7 where the increased power of (σ − 1) near to the true threshold for the second linear distribution means that the strength of the induced flow always remains smaller than that for the first linear distribution and for the uniform distribution.

Quadratic distribution
Given that the initial rise in the flow immediately after σ rises above unity is proportional to (σ − 1) 3 in the case of a uniform distribution and is proportional to (σ − 1) 4 for the second linear distribution, we have chosen to consider the following quadratic probability density function, This has been chosen because the powers of (σ − 1) quoted above seem to be associated with the behaviour of the probability density function at γ = 1. For the uniform distribution, Φ(1) is nonzero, while for the second linear distribution it is zero, but the slope is nonzero. The quadratic distribution has been chosen to have a zero value and a zero slope but a nonzero second derivative. We obtain the following, where K = 1 20 Clearly, we have obtained an initial velocity which is proportional to (σ −1) 5 as σ rises above unity. This further increased power of (σ − 1) is evident in Fig. 7. The pseudo-threshold is now which represents another increase over the two linear distribution cases; see Fig. 7b, d.

Seminormal distribution
We select to use the seminormal distribution which has a mean of unity (and hence a standard deviation which is also unity), and therefore the probability density function is Upon following the same ideas as above, we obtain the following expressions for the scaled velocity, permeability and porosity: The variation of the scaled velocity with the dimensionless pressure gradient, σ , is shown in Fig. 7. We have also found the following asymptotic expressions for F(σ ) for both large and small values of σ : For this distribution, the pseudo-threshold is This is clearly less than the mean of the probability density function, but it reflects the fact that flow is initiated, albeit weakly, at values of σ which are less than unity. The fact that this is a distribution with a mean of unity tells us that a proportion of the channels have widths which are greater than h and therefore flow may arise when σ < 1. For very small values of σ , Eq. (46) shows that the strength of the flow is superexponentially small. When σ = 1, we have F(σ ) = 0.360539, at which point all the other distributions which we have considered have just reached their common threshold values of σ .

Conclusions
It has been found that there is a very wide variety of relationships between the Darcy velocity and the applied pressure gradient when a porous medium is composed of distributions of channels or pipes. The chief difference between most of the different configurations of channel which are considered here is the manner in which flow begins once the scaled pressure gradient, σ , increases past unity. For a finite number of channels comprising a periodic configuration, the initial growth in the strength of the flow is quadratic in (σ − 1).
Once continuous distributions are encountered, the threshold becomes weaker, and the power of (σ − 1) which may be expected then depends on the exact nature of the probability density function, Φ(γ ), as γ → 1 − assuming that Φ(γ ) = 0 when γ > 1. We have seen two cases which are cubic, and have also obtained both quartic and quintic variations in (σ − 1). No doubt other distributions could be devised which would yield other types of variation. Specifically, if we restrict ourselves to a distribution, Φ(γ ), where Φ = 0 when γ > 1, then if the limit as γ → 1 − of the nth derivative of Φ(γ ) is nonzero, but that all other lower derivatives are zero, then the initial rise in F(σ ) at threshold will be proportional to (σ − 1) n+2 . Figure 7 also indicates that the pseudo-threshold displays a tendency to increase as the power of (σ − 1) increases. For the seminormal distribution, flow begins immediately once σ increases from zero, but its strength is super-exponentially small at first. Therefore, it is clear that the writing down of a definitive Darcy-Bingham model cannot be achieved in general, at least for a bundle of channels or tubes, but that the appropriate model will be specific to the type of porous medium being considered. This conclusion confirms the view of Coussot (2014) stated earlier. While upscaling techniques are often used very successfully in many other situations, the presence of yield surfaces at microscopic lengthscales, and the presence of stagnation in some pores and not others, would therefore seem to militate against the use of traditional upscaling to obtain what might be termed a universal Darcy-Bingham law.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
The aim of this "Appendix" section is to present a list of formulae for those cases where flow takes place in microchannels of circular cross section, namely Hagen-Poiseuille flow. We will assume that these channels have radius, h, and that H is the period for the underlying pattern of the different channels. Thus, for a single channel in a repeating square (H × H ) domain, the porosity is φ = π h 2 /H 2 , for reference. The formula for the fluid flux in a single channel which is analogous to that in Eq. (6) is given by This formula for g(σ ) is the classical Buckingham-Reiner formula for the flow of a Bingham fluid relative to that of a Newtonian fluid. Thus, Eqs. (9) and (8) may be replaced by, where K = φh 2 8 and φ = π h 2 H 2 .
The pseudo-threshold is at For the quadratic distribution the following may be found.