Effective rheology of two-phase flow in a capillary fiber bundle model

We investigate the effective rheology of two-phase flow in a bundle of parallel capillary tubes carrying two immiscible fluids under an external pressure drop. The diameter of each tube varies along its length and the corresponding capillary threshold pressures are considered to be distributed randomly according to a uniform probability distribution. We demonstrate through analytical calculations that a transition from a linear Darcy regime to a non-linear behavior occurs while decreasing the pressure drop $\Delta P$, where the total flow rate $\langle Q \rangle$ varies with $\Delta P$ with an exponent $2$. This exponent for the non-linear regime changes when a lower cut-off $P_m$ is introduced in the threshold distribution. We demonstrate analytically that, in the limit where $\Delta P$ approaches $P_m$, the flow rate scales as $\langle Q \rangle \sim (|\Delta P|-P_m)^{3/2}$. We have also provided some numerical results in support to our analytical findings.

Understanding the hydrodynamic properties of simultaneous flow of two or more immiscible fluids is essential due its relevance to a wide variety of different systems in industrial, geophysical and medical sectors [1,2]. Different applications, such as bubble generation in microfluidics, blood flow in capillary vessels, catalyst supports used in the automotive industry, transport in fuel cells, oil recovery, ground water management and CO 2 sequestration, deal with the flow of bubble trains in different types of systems, ranging from single capillaries to more complex porous media. The underlying physical mechanisms in multiphase flow are controlled by a number of factors, such as the capillary forces at the interfaces, viscosity contrast between the fluids, wettability and geometry of the system, which make the flow properties different from single phase flow. When one immiscible fluid invades a porous medium filled with another fluid, different types of transient flow patterns, namely viscous fingering [3,4], stable displacement [5], and capillary fingering [6] are observed while tuning the physical parameters [7]. These transient flow patterns were modeled by invasion percolation [8] and diffusion limited aggregation (DLA) models [9]. When steady state sets in after the initial instabilities, the flow properties in are characterized by relations between the global quantities, such as flow rate, pressure drop and fluid saturation [10,11]. It has been observed theoretically and experimentally that, in the regime where capillary forces compete with the viscous forces, the two-phase flow rate of Newtonian fluids in the steady state no longer obeys the linear Darcy law [12,13] but varies as a power law with the applied pressure drop [14][15][16][17]. Tallakstad et al. [14,15] experimentally measured the exponent of the power law to be close to two (= 1/0.54) in a two-dimensional system and followed this observation up with arguments why the exponent should be two. Rassi et al. [16] found a value for the exponent varying between 2.2 (= 1/0.45) and 3.0 (= 1/0.33) in a three-dimensional system. Sinha et al. [17] considered a similar system to that which had been studied by Rassi et al. finding an exponent 2.17±0.24 (= 1/(0.46±0.05)). The reason behind the discrepancy between the results of Rassi et al. and those of Sinha et al. is the possibility of a non-zero threshold pressure that observed in the later study, under which there would be no flow, which was assumed to be zero in the former study. The reciprocals in the brackets are provided in order to compare the exponent values reported in the literature [14][15][16][17] with those we present here in this article, as we express our results as Q as a power law in P, whereas in the cited papers P was expressed as a power law in Q .
This power law behavior is in contrast to the assumption of linearity in the relation between flow rate and pressure drop that is generally assumed in the relative permeability approach dominating reservoir simulations [18].
For a single capillary tube with varying diameter, Sinha et al. [19] showed that the average volumetric flow rate q in the steady state has a non-linear square-root type relationship with the pressure drop P as q ∼ P 2 − P 2 c . This was shown analytically by integrating the instantaneous linear two-phase flow equation over the whole capillary tube. Here P c is the threshold pressure difference below which there is no flow. It appears due to the capillary barriers at the interfaces at the narrow pore throats. Extending this non-linear relationship to a network of disordered pores, the relationship between the steady-state flow rate and an excess pressure drop leads to a quadratic relationship in the capillary dominated regime [20]. The quadratic relationship for the pore network, both in two and in three dimensions, was obtained analytically by mean-field calculations and numerically with pore network modeling [17,20].
While increasing the pressure drop, the capillary forces become negligible compared to the viscous forces. This leads to a crossover from the non-linear regime to a linear Darcy regime for both the single capillary tube and for the pore network. Such non-linear quadratic relationship at low flow rate and a crossover to a linear regime at high flow rate was also observed in case of the single-phase flow of Bingham viscoplastic fluid in porous media [21,22]. A Bingham fluid is a yield threshold fluid which behaves like a solid below the threshold and flows like a Newtonian fluid above it. The origin of the quadratic regime for the Bingham fluid flowing in a porous media can be understood intuitively in this way: the flow starts when one connected channel appears in the system just above a threshold pressure and the flow rate varies linearly with the excess pressure drop; while increasing the applied pressure drop further, more number of connected flow channels start to appear enhancing the overall flow rate more rapidly than the applied pressure drop leading to the quadratic relationship. Finally, when all possible flow paths become active, the flow become Newtonian following the linear Darcy law. Note that, in general, the rheology of the Bingham fluid is linear above the yield threshold. It is the disorder in the yield thresholds due to the porous medium that creates the quadratic regime.
The argument presented by Tallakstad et al. [14,15] focused on the successive opening of fluid channels when the pressure drop across the system was increased. When | P| is small, the flow will occur along isolated channels. The volumetric flow rate in such a channel will be proportional to | P|/L. Between the channels there will be fluid clusters held in place by capillary forces, say of the order p t . There is a pressure gradient | P|/L in the flow direction. A given cluster of length l will be stuck if p t > l | P|/L. The largest stuck cluster will then have a size l ,max = Lp t /| P|. If we now assume that this length, l ,max is same as the distance between the channels where there is flow, l ⊥ , then the total flow rate must be equal to the number of channels, which is proportional to 1/l ⊥ , multiplied by the flow rate in each channel. Hence, we have Q ∝ (1/l ⊥ ) | P| ∝ | P| 2 . Though this argument provides the same behavior as the one based on the mean field calculation [20] for two-dimensional networks, a difference appears in three dimensions. When following the same arguments, it leads the flow rate to vary with the pressure drop with third power as long as the isolated channels remain one-dimensional strings rather than two-dimensional sheets in three dimensions.
We present in this article a capillary fiber bundle model [23,24], which is a system of N parallel capillary tubes, disconnected from each other, each carrying an independent bubble trail of two immiscible fluids under an external pressure drop. In a porous medium, a typical pore consists of two wide pore bodies at the ends and a narrow pore throat in the middle. When an interface moves along the pore, the capillary pressure associated with the interface becomes position dependent due to the change in the radius of curvature. This introduces an overall threshold pressure that depends on the position of all the interfaces [19]. One can simplify the shape of the pore by a sinusoidal type and a long capillary tube with varying radius can be seen as a series of many pores. In the capillary bundle model, the diameter of each tube varies along the axis identically and the disorder in the threshold appear due to the different interface positions in different tubes. This model is essentially the only model for immiscible two-phase flow which is analytically tractable. We calculate the total average flow rate as a function of the applied pressure drop and study the effect of disorder in the threshold distribution. We point out that, here we do not address the question of the relation between the fluid distributions in the capillaries and the respective threshold distributions. Our aim with this model is to investigate how the range of the disorder in the threshold distribution controls the effective flow properties. This provides an insight into the non-linearities in steady-state two-phase flow. We will see that the exponent for the non-linear regime depends on the lower cut-off of the threshold distribution as well as on the behavior of the distribution near the cut-off. The possibility to study analytically for this model how the competition between viscous and capillary forces renders the Darcy relation non-linear, is a new and useful discovery.
The capillary fiber bundle model is a hydrodynamic analog of the fiber bundle model used in fracture mechanics to study mechanical failure under stress [25]. The fiber bundle model is an ideal example of a disordered system in statistical mechanics that is driven by threshold activated dynamics. It is a simple, yet very rich model to understand failure events in mechanical systems. In its simplest form it is analytically tractable. In more complex versions of the model, analytical calculations go hand in hand with numerical simulations. Figure 1 illustrates a bundle with N = 5 parallel capillaries. Each capillary tube has a length L and an average inner area a. For each capillary, the diameter varies along the long axis identically. Each capillary is filled with a bubble train of wetting and nonwetting fluids. Due to the varying diameter, the capillary forces at the interfaces vary as the bubble train moves along the tubes. We assume that the wetting fluid does not form films along the pore walls so that the fluids do not pass each other. The lengths of the wetting and non-wetting fluids in each tube is L w and L n , respectively such that the volume of the wetting fluid in each tube is L w a and the volume of the non-wetting fluid is L n a, where a is the average cross-sectional area of the capillary tubes. Hence, the saturations are given by S w = L w /L and S n = L n /L for each capillary tube. The cross-sectional pore area of the capillary fiber bundle is Though each tube contains the same amount of each fluid, it has its own division of the fluids into bubbles. We average over the ensemble of capillary tubes in the bundle by averaging over the fluids in each tube that pass at a given instance through an imaginary cut as shown in the figure. We will obtain the same averages if we consider a single capillary tube, averaging over a time interval the fluid passing the imaginary cut [19,26]. This shows that the model is ergodic. The volumetric flow rate in a capillary tube is given by [19]  where P is the pressure drop across the capillary tube, P c the sum of all the capillary forces along the capillary tube due to the interfaces and is the effective viscosity.
(| P| − P c ) is the Heaviside function which is zero for negative arguments and one for positive arguments.
Sinha et al. [19] showed that the time average when the pressure difference across the tube is kept fixed is given by where the function sgn( P) is the sign of the argument. Suppose now that the thresholds P c are distributed uniformly between zero and a maximum value P M . The cumulative threshold probability is then We have N capillary tubes. Using order statistics, we may order the N averaged threshold values, where 1 ≤ k ≤ N. Hence, The average volumetric flow rate through the capillary fiber bundle for | P| > 0 is then We assume the limit N → ∞ turning the sum into an integral, This integral is doable and we find when | P| ≤ P M and when | P| > P M . In the limit | P| ≫ P M , Equation (11) gives Hence, the Darcy relation for a tube is recovered. We see that this picture is consistent with that central to the arguments of Tallakstad et al. [14,15] leading to the quadratic dependence of Q on P. From Equation (7) we deduce that a number k c of the capillary tubes are active, where The typical distance between active capillary tubes in units of the distance between the tubes is then given by in accordance with the argument of Tallakstad et al. How stable is the square law Q ∝ | P| 2 ? That is, how much does it hinge on the choice of cumulative threshold probability (P c ). So far we have only considered the one given in Equation (5). Let us now generalize it to where α > 0. The average ordered threshold are then given by and when combined with the expression for Q , Equation (8) in the limit N → ∞, we find Since we are interested in the behavior for | P| → 0, we do this integral under the assumption that | P| < P M finding where the Ŵ function for real positive z is defined as, Ŵ(z) = ∞ −∞ t z−1 e −t dt. When α = 1, we recover Equation (10).
Equation (10) shows the behavior observed experimentally in References [14] and [15]. With Equation (18), we have just shown that Q /N ∼ | P| γ as | P| → 0, where γ depends on the threshold distribution, i.e., on α in Equation (15). Does this imply that there is no universality; that the experimentally observed behavior is due to the presence of a very specific threshold distribution?
As we now argue, there is universality. We note that the threshold distribution p(P c ) = d (P c )/dP c behaves as p(P c ) ∝ P α−1 c . Hence, if α > 1, the distribution vanishes as P c → 0, whereas it diverges for α < 1. Thus, the behavior of the distribution is vastly different for these two cases, and this causes γ to depend on α. However, for α = 1, the distribution reaches a constant, non-zero value for P c → 0. Any threshold distribution with this behavior for small P c , i.e., p(P c ) reaching a non-zero value and dp(P c )/dP c → 0 in the limit P c → 0 will give rise to the square power law seen in Equation (10). Such distributions are ubiquitous, and γ = 2 is universal over this class of distributions.
We now consider α = 1 again, but introduce a minimum threshold P m so that the cumulative threshold probability is given by Equation (6) yields in this case the ordered threshold sequence Equation ( Frontiers in Physics | www.frontiersin.org when we assume P m ≤ | P| ≤ P M . We find We find to lowest order in (| P| − P m ), that this expression behaves as as | P| → P m .
We now turn to numerical simulations and observe that the numerical results are in good a agreement with the analytical findings. The numerical simulations also allow us to explore the regions which are analytically challenging. Results are shown in Figure 2 for a bundle containing N = 10 5 capillary tubes and averaged over 10 4 configurations. In Figure 2A, we show the behavior of the volumetric flow rate Q as a function of increasing pressure drop P for uniform threshold distributions with P m = 0 and P m > 0, given by Equations (5, 19), respectively. The results show that, for each threshold distribution, the relationship is linear for high P obeying the Darcy law as predicted by Equation (12). For small pressure drops, Q follows a power law in P with an exponent 2 when there is no lower cutoff in the threshold distribution, i.e., P m = 0. This is predicted in Equation (10). When a lower cut-off is introduced in the threshold distribution (P m > 0), this exponent shifts from 2 to Frontiers in Physics | www.frontiersin.org 3/2 as predicted in Equation (23). These exponents in the nonlinear regime are not sensitive to the span of the distribution as shown in Figure 2A. An insight to a more generalized picture is presented in Figures 2B,C for a threshold distribution given by a generalization of Equation (15) with an introduction of a lower cut-off P m , With this distribution of thresholds, the exponent γ in the nonlinear region shows a continuous variation with α as γ = α + 1 for P m = 0. Such variation is given in Equation (18) and matches well with the numerical results as shown in Figure 2D.
In presence of a lower cut-off P m > 0, γ varies as (α+1/2) instead of (α + 1) irrespective of the position of the lower cut-off. An analytical treatment for a general α value with P m > 0 is rather challenging. Nevertheless, our numerical result matches with the analytical study (see Equation 23) in the limit α = 1. Equation (18) predicts an exponent γ = α + 1. A simple argument, related to that given by Roux and Herrmann [22], goes as follows: The number of active capillary tubes is proportional to (| P| − P m ) α . This behavior is observed in the insets in Figures 2B,C. The flow rate in an active capillary is proportional to (| P| − P m ) 1/2 . Hence, the total flow rate should be Q ∝ (| P| − P m ) α+1/2 . It is accidental that this argument works out for P m > 0 (Figure 2C), as it does not when P m = 0, where γ = α + 1. For the argument to function, the distribution of active capillaries and the flow rate in each capillary should be uncorrelated. It is not.
We find the same behavior with respect to the cut-off: An exponent 3/2 for the cumulative threshold probability where P m = 10 −β and P M = 10 β and β ranging from 0.5 to 1.5. The same goes for the cumulative threshold probability where we have set P m = 0.1 and P d = 1. In both of these cases, the probability density at P c = P m is finite.
We have presented an analytical study supported by numerical simulations of steady-state two-phase flow in a system of parallel capillary tubes. Considering a uniform distribution for the threshold pressures for the capillaries, we have calculated the average flow rate as a function of the applied pressure drop. When the thresholds are distributed according to a uniform distribution between zero and a maximum value-or more generally, the threshold distribution approaches a non-zero value in the limit of zero thresholds-we obtain a quadratic relationship between the flow rate and the applied pressure drop when the applied pressure drop is below the maximum threshold pressure, and the linear Darcy relationship for higher pressure drops. This crossover between a quadratic non-linear and linear flow regimes is in agreement with many existing results of two-phase flow in porous media which shows that this simple model can capture effective two-phase flow properties of more complex porous media. When a lower cut-off is introduced in the threshold distribution, the quadratic relationship changes, and the flow rate varies with an excess pressure drop with an exponent 3/2 as the pressure drop approaches to the lowest threshold pressure.
The difference between the capillary fiber bundle model and a porous medium is that in the latter, the fluids meet and mix at the nodes of the pore network. This is an essential mechanism that leads to the non-linear Darcy law is a power law with an exponent two as seen in the experiments, the numerical simulations and the mean-field calculations. However, it remains a mystery how the mixing at the nodes leads to this universality.

DATA AVAILABILITY
The datasets for this manuscript are not publicly available because all data sets are included in manuscript. Requests to access the datasets should be directed to subhadeep.roy@ntnu.no.

AUTHOR CONTRIBUTIONS
AH developed the theory and performed the analytical calculations of the manuscript. AH and SS wrote the first draft of the manuscript. SR did the numerical simulations. All the authors contributed in developing the theory and writing the manuscript to its final form.