Original Research ARTICLE
Effective Rheology of Two-Phase Flow in a Capillary Fiber Bundle Model
- 1PoreLab, Department of Physics, Norwegian University of Science and Technology (NTNU), Trondheim, Norway
- 2Beijing Computational Science Research Center, Beijing, China
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 the tubes vary along the length which introduce capillary threshold pressures. We demonstrate through analytical calculations that a transition from a linear Darcy to a non-linear behavior occurs while decreasing the pressure drop ΔP, where the total flow rate 〈Q〉 varies with ΔP with an exponent 2 as 〈Q〉~ΔP2 for uniform threshold distribution. The exponent changes when a lower cut-off Pm is introduced in the threshold distribution and in the limit where ΔP approaches Pm, the flow rate scales as . While considering threshold distribution with a power α, we find that the exponent γ for the non-linear regime vary as γ = α + 1 for Pm = 0 and γ = α + 1/2 for Pm > 0. We provide numerical results in support of 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 CO2 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 , and capillary fingering  are observed while tuning the physical parameters . These transient flow patterns were modeled by invasion percolation  and diffusion limited aggregation (DLA) models . 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–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. 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. 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–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 .
For a single capillary tube with varying diameter, Sinha et al. showed that the average volumetric flow rate in the steady state has a non-linear square-root type relationship with the pressure drop ΔP as . This was shown analytically by integrating the instantaneous linear two-phase flow equation over the whole capillary tube. Here Pc 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 . 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 pt. There is a pressure gradient |ΔP|/L in the flow direction. A given cluster of length l∥ will be stuck if pt > l∥|ΔP|/L. The largest stuck cluster will then have a size l∥, max = Lpt/|Δ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 . Though this argument provides the same behavior as the one based on the mean field calculation  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 . 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 . 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 non-wetting 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 Lw and Ln, respectively such that the volume of the wetting fluid in each tube is Lwa and the volume of the non-wetting fluid is Lna, where a is the average cross-sectional area of the capillary tubes. Hence, the saturations are given by Sw = Lw/L and Sn = Ln/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.
Figure 1. The capillary tube model. There are N = 5 capillaries, each filled with a bubble train of wetting (white) and non-wetting (black) fluids moving in the direction of the arrow. The diameter along each tube vary so that the capillary force from each interface vary with its position. The variation in the diameters are not illustrated in the figure. The average diameters are the same for all tubes. The broken line illustrates an imaginary cut through the capillary fiber bundle.
The volumetric flow rate in a capillary tube is given by 
where ΔP is the pressure drop across the capillary tube, Pc the sum of all the capillary forces along the capillary tube due to the interfaces and
is the effective viscosity. Θ(|ΔP| − Pc) is the Heaviside function which is zero for negative arguments and one for positive arguments.
Sinha et al. 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 Pc are distributed uniformly between zero and a maximum value PM. 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| ≤ PM and
when |ΔP| > PM. In the limit |ΔP| ≫ PM, 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 kc 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 Π(Pc). 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| < PM finding
where the Γ function for real positive z is defined as, . When α = 1, we recover Equation (10).
Equation (10) shows the behavior observed experimentally in References  and . 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(Pc) = dΠ(Pc)/dPc behaves as . Hence, if α > 1, the distribution vanishes as Pc → 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 Pc → 0. Any threshold distribution with this behavior for small Pc, i.e., p(Pc) reaching a non-zero value and dp(Pc)/dPc → 0 in the limit Pc → 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 Pm so that the cumulative threshold probability is given by
Equation (6) yields in this case the ordered threshold sequence
Equation (8) now becomes in the limit N → ∞
when we assume Pm ≤ |ΔP| ≤ PM. We find
We find to lowest order in (|ΔP| − Pm), that this expression behaves as
as |ΔP| → Pm.
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 = 105 capillary tubes and averaged over 104 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 Pm = 0 and Pm > 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 cut-off in the threshold distribution, i.e., Pm = 0. This is predicted in Equation (10). When a lower cut-off is introduced in the threshold distribution (Pm > 0), this exponent shifts from 2 to 3/2 as predicted in Equation (23). These exponents in the non-linear 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 Pm,
With this distribution of thresholds, the exponent γ in the non-linear region shows a continuous variation with α as γ = α + 1 for Pm = 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 Pm > 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 Pm > 0 is rather challenging. Nevertheless, our numerical result matches with the analytical study (see Equation 23) in the limit α = 1.
Figure 2. Results from numerical simulations performed with N = 105 capillaries and averaged over 104 configurations. Variation of 〈Q〉/N as a function of pressure drop ΔP for different threshold distributions are shown in (A–C) where non-linear to linear transitions are observed while increasing the pressure drop. (A) Corresponds to uniform threshold distribution (Equation 5) where the power-law exponent γ for the non-linear regime has a value 2 without a lower cut-off (Pm = 0). With any non-zero lower cut-off (Pm > 0), the exponent shifts to 3/2 (Equation 23). Results for the threshold distribution with a power α (Equation 24) are shown in (B,C) for Pm = 0 and Pm > 0, respectively, where γ varies with α as γ = α + 1 for Pm = 0 and as γ = α + 1/2 for Pm > 0. These two relations are shown in (D) for the range of α which show two distinct straight lines for Pm = 0 and for Pm > 0. Here, the number of active capillaries (kc) vary with ΔP as as shown in the insets of (B,C).
Equation (18) predicts an exponent γ = α + 1. A simple argument, related to that given by Roux and Herrmann , goes as follows: The number of active capillary tubes is proportional to . This behavior is observed in the insets in Figures 2B,C. The flow rate in an active capillary is proportional to . Hence, the total flow rate should be . It is accidental that this argument works out for Pm > 0 (Figure 2C), as it does not when Pm = 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 and and β ranging from 0.5 to 1.5. The same goes for the cumulative threshold probability
where we have set Pm = 0.1 and Pd = 1. In both of these cases, the probability density at Pc = Pm 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.
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 firstname.lastname@example.org.
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.
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 thank Dick Bedeaux, Carl Fredrik Berg, Eirik G. Flekkøy, Signe Kjelstrup, Knut Jørgen Måløy, Per Arne Slotte, and Ole Torsæter for interesting discussions. This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. SS was supported by the National Natural Science Foundation of China under grant number 11750110430.
10. Valavanides M. Review of steady-state two-phase flow in porous media: independent variables, universal energy efficiency map, critical flow conditions, effective characterization of flow and pore network. Transp Porous Media. (2018) 123:45–99. doi: 10.1007/s11242-018-1026-1
11. Hansen A, Sinha S, Bedeaux D, Kjelstrup S, Gjennestad MA, Vassvik M. Relations between seepage velocities in immiscible, incompressible two-phase flow in porous media. Transp Porous Media. (2018) 125:565–87. doi: 10.1007/s11242-018-1139-6
12. Darcy H. Les Fontaines Publiques de la Ville de Dijon: Exposition et Application des Principes à Suivre et des Formules à Employer dans les Questions de Distribution d'eau. Paris: V. Dalamont (1856).
14. Tallakstad KT, Knudsen HA, Ramstad T, Løvoll G, Måløy KJ, Toussaint R, et al. Steady-state two-phase flow in porous media: statistics and transport properties. Phys Rev Lett. (2009) 102:074502. doi: 10.1103/PhysRevLett.102.074502
15. Tallakstad KT, Løvoll G, Knudsen HA, Ramstad T, Flekkøy EG, Måløy KJ. Steady-state, simultaneous two-phase flow in porous media: an experimental study. Phys Rev E. (2009) 80:036308. doi: 10.1103/PhysRevE.80.036308
16. Rassi EM, Codd SL, Seymour JD. Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow. New J Phys. (2011) 13:015007. doi: 10.1088/1367-2630/13/1/015007
17. Sinha S, Bender AT, Danczyk M, Keepseagle K, Prather CA, Bray JM, et al. Effective rheology of two-phase flow in three-dimensional porous media: experiment and simulation. Transp Porous Media. (2017) 119:77–94. doi: 10.1007/s11242-017-0874-4
21. Chevalier T, Talon L. Generalization of Darcy's law for Bingham fluids in porous media: from flow-field statistics to the flow-rate regimes. Phys Rev E. (2015) 91:023011. doi: 10.1103/PhysRevE.91.023011
Keywords: two-phase flow, capillary fiber bundle model, effective rheology, non-Darcy flow at low velocity, porous media
Citation: Roy S, Hansen A and Sinha S (2019) Effective Rheology of Two-Phase Flow in a Capillary Fiber Bundle Model. Front. Phys. 7:92. doi: 10.3389/fphy.2019.00092
Received: 23 February 2019; Accepted: 11 June 2019;
Published: 09 July 2019.
Edited by:Ferenc Kun, University of Debrecen, Hungary
Reviewed by:Allbens Picardi Faria Atman, Federal Center for Technological Education of Minas Gerais, Brazil
Muktish Acharyya, Presidency University, India
Copyright © 2019 Roy, Hansen and Sinha. 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) and the copyright owner(s) 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: Subhadeep Roy, email@example.com