Steady-State Two-Phase Flow of Compressible and Incompressible Fluids in a Capillary Tube of Varying Radius

Abstract

We study immiscible two-phase flow of a compressible and an incompressible fluid inside a capillary tube of varying radius under steady-state conditions. The incompressible fluid is Newtonian and the compressible fluid is an inviscid ideal gas. The surface tension associated with the interfaces between the two fluids introduces capillary forces that vary along the tube due to the variation in the tube radius. The interplay between effects due to the capillary forces and the compressibility results in a set of properties that are different from incompressible two-phase flow. As the fluids move towards the outlet, the bubbles of the compressible fluid grows in volume due to the decrease in pressure. The volumetric growth of the compressible bubbles makes the volumetric flow rate at the outlet higher than at the inlet. The growth is not only a function of the pressure drop across the tube, but also of the ambient pressure. Furthermore, the capillary forces create an effective threshold below which there is no flow. Above the threshold, the system shows a weak nonlinearity between the flow rates and the effective pressure drop, where the nonlinearity also depends on the absolute pressures across the tube.

1 Introduction

Hydrodynamic properties of the flow of multiple immiscible and incompressible fluids, otherwise known as two-phase flow (Bear 1988; Dullien 1992; Blunt 2017; Feder et al. 2022), are controlled by a number of different factors: fluid properties such as the viscosity contrast and surface tension between the fluids, driving parameters such as the applied pressure drop or the flow rate, and geometrical properties such as the size and shape of the space in which the fluids are flowing. The combined effects of these factors make two-phase flow different and more complex than single phase flow. The dimensionless parameters that play a key role to define the flow properties are the ratio between the viscous and capillary forces, referred to as the capillary number, and the ratio between the viscosities of the two fluids. Depending on the values of these parameters, the flow generates different types of fingering patterns (Chen and Wilkinson 1985; Lenormand and Zarcone 1985; Måløy et al. 1985; Løvoll et al. 2004; Zhao et al. 2019) or stable displacement fronts (Lenormand and Touboul 1988) during invasion processes where one fluid displaces another in a porous medium.

Displacement processes are transient. If one continues to inject after breakthrough, the flow enters a steady state characterized by a situation where the macroscopic flow properties fluctuate or remain constant around well-defined averages. A more general form of steady-state flow can be achieved by continuously injecting both fluids simultaneously. In this case, the dynamics at the pore scale might have fluid clusters breaking up and forming, while the macroscopic flow parameters still have well-defined averages.

Over the last decade, it has become clear that steady-state flow deviates from the linear Darcy relationship (Darcy 1856) between the total flow rate and pressure drop over a range of parameters. Rather, one finds a power law relationship between pressure drop and the volumetric flow rate (Tallakstad et al. 2009; Rassi et al. 2011) in that range. In terms of the capillary number, this range is intermediary, with linearity appearing both for lower and higher values (Sinha et al. 2017; Gao et al. 2020; Zhang et al. 2021). Theoretical work to understand the physics behind the nonlinearity has appeared in, e.g., Tallakstad et al. (2009); Sinha and Hansen (2012); Zhang et al. (2021), and computational studies have been performed using Lattice Boltzmann simulations (Yiotis et al. 2013) and dynamic pore network modeling (Sinha et al. 2021, 2017). It is now believed that a fundamental mechanism behind this nonlinearity is the capillary barriers at the pore throats, which create an effective yield threshold. When the viscous forces increase, they overcome the capillary barriers creating new flow paths. This increases the effective mobility and thus the nonlinear behavior appears (Roux and Herrmann 1987). The disorder in the pore-space properties, such as the pore-size distribution (Roy et al., 2021) and the wetting angle distribution (Fyhn et al., 2021), therefore play key roles in determining the value of the exponent relating the volumetric flow rate and the pressure drop in the nonlinear regime.

The majority of the analytical and numerical approaches mentioned above consider the two fluids to be incompressible, whereas many of the experiments and applications use air as one of the fluids. Air is strongly compressible, which can enhance the complex pore-scale mechanisms such as trapping and coalescence (Leverett 1941; Li and Yortsos 1994). Compressibility is relevant to a wide range of applications with liquid and gas transport in porous media, for example, $\textrm{CO}_2$ transport and storage (Reynolds and Krevor 2015; Abidoye et al. 2015; Iglauer et al. 2019) and the transport in fuel cells (Niblett et al. 2020). Another class of applications where the compressibility plays a key role are those involving phase transitions of the fluids such as boiling and condensation. There are industrial applications where such processes are of high importance, for example aerospace vehicle thermal protection (Huang et al. 2017), high power electronics cooling systems (Gedupudi et al. 2011; Li et al. 2012, 2020) and chemical reactors (Bremer and Sundmacher 2019). These applications utilize the high specific surface area of a porous medium with fluid flowing inside, which enhances the heat and mass transfer rates (Sapin et al. 2016; Sun et al. 2011). There are also natural processes such as drying of soil (Rossi and Nimmo 1994) where a liquid to gas transition takes place.

In this article, we present a study of two-phase flow of a mixture of compressible and incompressible fluids in a capillary tube with varying radius. We consider two fluids, one is an incompressible Newtonian fluid obeying Poiseuille flow in the steady state whereas the other is a compressible ideal gas, where the viscosity is assumed to be negligible. The fluids flow as a series of bubbles and droplets under a constant pressure drop along the tube.

In case of two-phase flow of two incompressible fluids in a corresponding capillary tube, it has been found that the volumetric flow rate (Q) depends on the square root of the pressure drop ($\Delta P$) along the tube minus a threshold pressure ($P_t$), that is, $Q\sim \sqrt{\Delta P - P_t}$ (Sinha et al. 2013). One primary goal of the present work is to determine how this constitutive equation changes when one of the two fluids is compressible.

A secondary goal of this work is to provide a basis for dynamic pore network modeling (Blunt 2001; Meakin and Tartakovsky 2009; Joekar-Niasar and Hassanizadeh 2012; Sinha et al. 2021) of compressible-incompressible fluid mixtures. This opens the possibility for incorporating thermodynamic effects in such models such as boiling. However, in order to explore the effect of compressibility on the rheological properties in general, we considered a higher range of pressure drops here, whereas some specific applications mentioned earlier in this section may need a different range. We also note that the other dominating computational model in this context, the Lattice Boltzmann model (Gunstensen et al. 1991; Ramstad et al. 2012), can only incorporate fluids that are weakly compressible (Qiu et al. 2017; Guo et al. 2020).

We describe in Sect. 2.1 the equations that govern the flow through the capillary tube. In Sect. 2.2, we introduce the boundary conditions used, i.e., how we inject alternate compressible and incompressible fluid into the tube. Sect. 2.3 describes how the governing equations are integrated in time.

Section 3 presents the results of our investigation. Section 3.1 defines what we mean by steady-state flow in the context of expanding bubbles. In Sect. 3.2, we investigate how the compressible bubbles grow as they advance along the tube, thus increasing the overall flow rate of the fluids. Section 3.3 presents the relation between volumetric flow rate and pressure drop at both the inlet and outlet.

We summarize our results in Sect. 4. Section 5 contains the description of the videos provided in the electronic supplementary material.

2 Methodology

The capillary tube considered in this work is filled with an incompressible and a compressible fluid, immiscible to each other, which flow through it. The fluids are separated by menisci associated with a surface tension. In order to introduce a variation in the capillary forces along the tube, we consider a periodic variation in the radius of the capillary tube along the flow direction x. The incompressible fluid is a viscous Newtonian liquid obeying Hagen-Poiseuille flow whereas the compressible fluid is an inviscid ideal gas. The flow occurs as a plug flow with a series of alternate bubbles and droplets of the two fluids as illustrated in Fig. 1. There is no fluid film along the tube walls and therefore no coalescence or snap off taking place inside the tube during the flow. We will refer the compressible and the incompressible fluid segments as bubbles and droplets, respectively.

Fig. 1

Illustration of the tube geometry and the indexed variables. The shaded fluid represents the non-wetting compressible gas and the white fluid represents the wetting incompressible liquid. There are $N=6$ bubbles here indicated by the numbers $i=1,\ldots ,6$. The indexed variables $P_i$, $V_i$ and $n_i$, respectively, correspond to the pressure, volume and moles of the ith bubble whereas $Q_i$ corresponds to the flow rate of the droplet between ith and $(i+1)$th bubbles.

2.1 Governing Equations

We assume that at a given time the system contains N compressible bubbles denoted by $i=1,2,\ldots ,N$ from left to right as shown in Fig. 1. The volume $V_i$ and the pressure $P_i$ of the ith bubble are connected through the ideal gas law,

$$\begin{aligned} \displaystyle P_iV_i = n_iRT \;, \end{aligned}$$

(1)

where $n_i$ is the number of moles of gas present inside the bubble, R is the ideal gas constant and T is the temperature. The volume of an incompressible droplet on the other hand will remain constant throughout the flow and the flow rate will depend on the pressures of the two compressible bubbles bordering it. The volumetric flow rate of the incompressible droplet between i and $i+1$ is denoted by $Q_i$, and follows the constitutive equation (Dullien 1992; Washburn 1921),

$$\begin{aligned} \displaystyle Q_i = \frac{A^2}{8\pi \mu (x_{i+1}^l-x_{i}^r)}\left[ P_i - P_\textrm{c}(x_i^r) - P_{i+1} + P_\textrm{c}(x_{i+1}^l)\right] \;, \end{aligned}$$

(2)

where $\mu$ is the viscosity of the incompressible fluid and $P_c(x)$ is the capillary pressure at x. Here we assumed that the variation in the tube radius only affects the capillary pressure $P_c$ along the tube and therefore the area A in the above equation is considered to be the average cross-sectional area of the tube. This is an approximation that is commonly used in dynamic pore network models (Sinha et al. 2021). Furthermore, the bubbles are assumed to be smaller in size compared to the period of the tube so that the flow of the incompressible bubbles in the slowly-changing area can be considered locally as a Poiseuille flow (Panton 2013). The volume of a compressible bubble is therefore given by, $V_i=A(x_i^r-x_i^l)$ where $x_i^l$ and $x_i^r$ are the positions of the left and right menisci of the ith bubble, respectively.

Here we consider the incompressible fluid to be more wetting with respect to the pore walls than the compressible fluid, thus determining the sign of $P_c$ in Equation 2. We model $P_c$ by using the Young-Laplace equation (Dullien 1992),

$$\begin{aligned} \displaystyle P_\mathrm{{c}}(x) = \frac{2\gamma }{r(x)} \;, \end{aligned}$$

(3)

where r(x) is the radius of the tube at x. Here $\gamma = \sigma \cos (\theta )$ where $\sigma$ is the surface tension between the fluids and $\theta$ is the wetting angle of the fluid with respect to the tube wall. The variation in the radius of the tube shown in Fig. 1 is modeled by

$$\begin{aligned} \displaystyle r(x) = \frac{1}{2}\left[ w + 2a\cos \left( \frac{2h\pi x}{L}\right) \right] \end{aligned}$$

(4)

where L is the tube length, w is the average radius, a is the amplitude of oscillation and h is the number of periods.

2.2 Boundary Conditions

The system is driven by a constant pressure drop $\Delta P = P_0-P_L$ where $P_0$ and $P_L$ are the pressures at the inlet ($x=0$) and outlet ($x=L$), respectively. The two fluids are injected alternatively at the inlet. Depending on the fluid that is being injected and the fluid that is leaving the tube, there will be different configurations as illustrated in Fig. 2. When a bubble is entering at the inlet [Figure 2(a)] or leaving at the outlet [Figure 2(c)], the pressure in that bubble is given by $P_0$ or $P_L$, respectively. This is because the compressible fluid has no viscosity and thus the pressure inside a bubble is uniform. The pressures inside all other bubbles are calculated using Equation 1. When a droplet is entering at the inlet [Figure 2(b) and (c)] or leaving at the outlet [Figure 2(a) and (b)], the respective flow rates $Q_0$ and $Q_N$ are given by,

$$\begin{aligned} \begin{aligned} \displaystyle Q_0&= \frac{A^2}{8\pi \mu x_1^l}\left[ P_0 - P_1 + P_\textrm{c}(x_1^l)\right] \; \textrm{and}\\ Q_N&= \frac{A^2}{8\pi \mu (L-x_N^r)}\left[ P_N - P_\textrm{c}(x_N^r) - P_L\right] \;, \end{aligned} \end{aligned}$$

(5)

whereas the flow rates of the remaining droplets are calculated using Equation 2.

Fig. 2

Illustration of different configurations where bubbles and droplets are colored as gray and white, respectively. In a, a bubble is entering at the inlet and therefore $P_1 = P_0$ there. In c, a bubble is leaving at the outlet, therefore $P_N=P_L$ there. A droplet is entering at the inlet in (b) and (c), and leaving at the outlet in a and b. The flow rates of such droplets are calculated using Equation 5.

The simulation is started with the tube completely filled with the incompressible fluid. The two fluids are then injected alternately at the inlet using small time steps. Whenever the injection is switched to a different fluid, a new menisci is created and the injection is continued for that fluid until the bubble or the droplet being injected has reached a given length, $b_\textrm{C}$ or $b_\textrm{I}$, respectively. For each new bubble or droplet, a new value for $b_\textrm{C}$ or $b_\textrm{I}$ is determined using the following scheme:

$$\begin{aligned} \displaystyle b_\textrm{C} = b_\textrm{min} + kF_\textrm{C}b_\textrm{max} \quad \; \textrm{and} \quad b_\textrm{I} = b_\textrm{min} + kF_\textrm{I}b_\textrm{max} \;, \end{aligned}$$

(6)

where k is chosen from a uniform distribution of random numbers between 0 and 1. $F_\textrm{C}$ and $F_\textrm{I}$ are the tentative values of the fractional flows for the bubbles and droplets, respectively. The two parameters $b_\textrm{max}$ and $b_\textrm{min}$ set the smallest and largest allowed sizes of any bubble or droplet. We consider here $b_\textrm{min}=L/10^4$ and $b_\textrm{max}=L/50$. The parameters $b_\textrm{C}$ and $b_\textrm{I}$ decide the initial sizes of the bubbles and droplets just after they detach from the inlet. For the compressible fluid, this determines the number of moles $n_i$ inside a bubble,

$$\begin{aligned} \displaystyle n_i = \frac{Ab_\textrm{C}P_0}{RT} \;, \end{aligned}$$

(7)

which remains constant for that bubble throughout the flow after it gets detached from the inlet.

2.3 Updating the Menisci Positions

At any time, the two menisci bordering a droplet inside the tube move with the same velocities. The velocities of the menisci are calculated from the velocities $v_i$ of the droplets using Equations 2 and 5,

$$\begin{aligned} \displaystyle \frac{{\hbox {d}x_{i}^r}}{\hbox {d}t} = \frac{{\hbox {d}x_{i+1}^l}}{\hbox {d}t} = v_i = \frac{Q_i}{A} \;. \end{aligned}$$

(8)

We solve these ordinary differential equations using an explicit Euler scheme, thus updating positions of all menisci by choosing a small time step $\Delta t$.

Depending on the positions of the menisci and the corresponding capillary pressures, the bubbles may compress or expand. If a bubble compresses at any time step, it means the left and right interfaces of that bubble approach each other. This necessitates the choice of time step $\Delta t$ to be sufficiently small, as otherwise, the two menisci around that bubble will collapse after the time step. We deal with this situation in the following way. First we calculate a time $\Delta t_1$ that is needed to pass one pore-volume of incompressible fluid through the tube,

$$\begin{aligned} \displaystyle \Delta t_1 = \frac{8\pi \mu L^2}{A(P_0-P_L)} \;. \end{aligned}$$

(9)

Next, we check for every bubble i if $(v_{i-1}-v_i)>0$, that is, whether the two menisci bordering the bubble are approaching each other in that time step. If this criterion is found to be true for any of the bubbles j, we measure the time it will take for the two menisci to collapse,

$$\begin{aligned} \displaystyle \Delta t_2^j = \frac{x_j^r-x_j^l}{v_{j-1}-v_j} \;. \end{aligned}$$

(10)

After calculating $\Delta t_1$ and $\Delta t_2^j$, we determine a time $\Delta t$ for that step from,

$$\begin{aligned} \displaystyle \Delta t = \textrm{min}(a^*\Delta t_1, b^*\Delta t_2^j) \;, \end{aligned}$$

(11)

which means that if there is a possibility for a bubble to collapse during the time step, we chose $\Delta t$ from the minimum of $a^*\Delta t_1$ and all of $b^*\Delta t_2^j$. If there is no possibility of collapse, we use $\Delta t$ equal to $a^*\Delta t_1$. For the simulations presented in this paper, we set $a^*=10^{-8}$ and $b^*=10^{-6}$.

3 Results and Discussions

We perform steady-state simulations considering a tube of length $L=100\,\textrm{cm}$ with $w=1\,\textrm{cm}$, $a=0.25\,\textrm{cm}$ and $h=30$ (Equation 4). The viscosity of the incompressible fluid is $\mu =0.001\,\mathrm{Pa.s}$, the ideal gas constant is $R=8.31\,\mathrm{J/(mol.K)}$ and the temperature is kept fixed throughout the simulation at $T=293\,\textrm{K}$. We fix $F_c=0.4$ (Equation 6) which sets the volumetric fractional flow of the compressible fluid at the inlet around that value. We perform simulations varying the pressure drops ($\Delta P = P_0-P_L$) as well as the absolute outlet pressure with different values of the surface tension, $\gamma$.

Fig. 3

Total volumetric flow rate $q_\textrm{T}^\textrm{i}$ at the inlet as a function of the injected pore volume $V_p$ for the outlet pressures (a) $P_L=1\,\textrm{kPa}$ and (b) $P_L=100\,\textrm{kPa}$. Here the surface tension $\gamma =0.09\,\mathrm{N/m}$. The steady-state values of the flow rates are measured by taking averages in the range of 20 to 40 pore volumes as indicated by the dashed lines. Here we only show the data sets corresponding to the pressure drops $\Delta P=1,2,4$ and $8\,\textrm{kPa}$ in order to keep the clarity, however all the data sets show the similar trend of reaching the steady state.

3.1 Steady-State Flow

The steady state is defined by the volumetric flow rates of the fluids fluctuating around a stable average. Due to the expansion of the compressible fluid, which we will discuss in a moment, the volumetric flow rate of the fluids changes as the fluids flow towards the outlet. We define the quantities $Q_\textrm{T}^\textrm{i}$, $Q_\textrm{C}^\textrm{i}$, $Q_\textrm{I}^\textrm{i}$ as the average steady-state flow rates for the total, compressible and incompressible fluids at the inlet and $Q_\textrm{T}^\textrm{o}$, $Q_\textrm{C}^\textrm{o}$, $Q_\textrm{I}^\textrm{o}$ as those at the outlet. The inlet and outlet flow rates are measured by tracking the displacements of the first meniscus nearest to the inlet and the last meniscus near the outlet, which are either the left or the right meniscus of the first ($i=1$) and the last ($i=N$) bubbles. The instantaneous flow rates of the bubbles and droplets are measured as $q_\textrm{C}^\textrm{i}=A\sum \Delta x_1^r/\sum \Delta t$ for $x_1^l=0$, $q_\textrm{I}^\textrm{i}=A\sum \Delta x_1^l/\sum \Delta t$ for $x_1^l>0$ and $q_\textrm{C}^\textrm{o}=A\sum \Delta x_N^l/\sum \Delta t$ for $x_N^l=L$, $q_\textrm{I}^\textrm{o}=A\sum \Delta x_N^r/\sum \Delta t$ for $x_N^R<L$. This measurement is performed after every 0.05 pore-volumes of fluid are injected and the sum is therefore over the time steps in between. The total flow rates are therefore given by, $q_\textrm{T}^\textrm{i,o} = q_\textrm{C}^\textrm{i,o}+q_\textrm{I}^\textrm{i,o}$. This provides the measurement of the injected and outlet flow rates as a function of the injected pore volumes or of the time. In Fig. 3, we plot $q_\textrm{T}^\textrm{i}$ as a function of the pore-volumes ($V_p$) injected for (a) $P_L=1\,\textrm{kPa}$ and (b) $P_L=100\,\textrm{kPa}$. The pore-volume $V_p$ is defined as the ratio between the total volume of the inject fluids and the volume of the total pore space of the tube, which provides an estimate of how many times the pore space was flushed with the fluids. All the plots show that the total flow rate $q_\textrm{T}^\textrm{i}$ increases with time at the beginning of the flow. This increase in $q_\textrm{T}^\textrm{i}$ is due to the decrease in the effective viscosity of the system caused by the injection of inviscid compressible gas into the tube filled with viscous incompressible fluid. After the injection of a few pore volumes, $q_\textrm{T}^\textrm{i}$ fluctuates around a constant average ($Q_\textrm{T}^\textrm{i}$) shown by the horizontal dashed lines which defines the steady state. Notice that, the averages are the same for the same values of $\Delta P$ for the two different outlet pressures $P_L$, however the fluctuations are different. This we will see more in the following, that the outlet pressure $P_L$ plays a significant role in the flow properties in addition to the pressure drop $\Delta P$. We run our simulations for 40 pore volumes of fluid where the steady-state averages are taken after 20 pore volumes injected to ensure that a steady state has been reached.

3.2 Bubble Growth

As a compressible bubble moves along the tube, the volume of the bubble increases due to the decrease in the pressure towards the outlet (Vazquez et al. 2010). The bubble can also grow due to other mechanisms, such as the increase in temperature or a phase transition between liquid and gas phases (Welch 1998; Kenning et al. 2006), but these phenomena are not studied here. A simulation with a single bubble inside a short tube is shown in the supplementary material which illustrates that the bubble increases in size as it flows towards the outlet. To understand how this growth depends on different flow parameters in the steady state, we define the growth function $G_\textrm{C}(x)$ by,

$$\begin{aligned} \displaystyle G_\textrm{C}(x) = \frac{V(x) - V_\textrm{0}}{V_0}\;, \end{aligned}$$

(12)

where $V_0$ and V(x) are the volume of a given bubble initially after detaching from the inlet and when its center is at x. We measure $G_\textrm{C}$ by including all the bubbles that are not attached to the inlet or outlet and calculate the time average value of $(V(x)-V_0)/V_0$ in the investigated time interval, where x is the center of the bubble.

Fig. 4

Plot of the bubble growth $G_\textrm{C}(x)$ in the steady state as a function of the scaled position x/L inside the tube for zero surface tension, $\gamma = 0$. The two plots show the results for the same set of pressure drops $\Delta P$ with different outlet pressures $P_L$.

Fig. 5

Variation of the pressure P(x) [kPa] inside a compressible bubble along the tube during steady state flow. P(x) shows a linear behavior for different values of $\Delta P$ and $P_L$.

Figure 4 shows the variation of $G_\textrm{C}(x)$ along the tube for two different outlet pressures, $P_L=1\,\textrm{kPa}$ and $100\,\textrm{kPa}$ where we plot the results for the same set of pressure drops $\Delta P$. These results are with zero surface tension, $\gamma =0$. There are a few details to note here. First, the plots show that $G_\textrm{C}(x)$ increases with an increase in $\Delta P$. In addition, $G_\textrm{C}(x)$ also depends on the absolute pressures at the inlet and outlet, since we can see that the curves are nonlinear functions of x for $P_L=1\,\textrm{kPa}$, whereas for $P_L=100\,\textrm{kPa}$, they show linear behavior. Furthermore, $G_\textrm{C}(x)$ approaches $\Delta P/P_L$ at $x=L$ for all the data sets.

To explain the dependency of $G_\textrm{C}(x)$ on $\Delta P$ and $P_L$, we recall Equation 1 and rewrite Equation 12 as,

$$\begin{aligned} \displaystyle G_\textrm{C}(x) = \frac{P_0-P(x)}{P(x)} \;, \end{aligned}$$

(13)

where P(x) is the pressure inside a bubble at x. For $x=L$, $P(x) = P_L$ and therefore $G_\textrm{C}(L)=\Delta P/P_L$ as observed. In Fig. 5, we plot P(x), averaged over different time steps in the steady state, for the two outlet pressures, $P_L=1\,\textrm{kPa}$ and $100\,\textrm{kPa}$. Both of the plots show linear variation along x with the slope $-\Delta P$. We therefore have $P(x)=-x\Delta P/L+P_L+\Delta P$ and thus,

$$\begin{aligned} \displaystyle \frac{G_\textrm{C}(x)}{n_P} = \frac{x/L}{1+n_P(1-x/L)} \;, \end{aligned}$$

(14)

where $n_P=\Delta P/P_L$. This leads to

$$\begin{aligned} \displaystyle \frac{G_\textrm{C}(x)}{n_P} \sim {\left\{ \begin{array}{ll} \displaystyle \frac{1}{n_P}\left( \frac{1}{1-x/L}-1\right) &{} \textrm{for} \quad n_P \gg 1 \;, \\ x/L &{} \textrm{for} \quad n_P \ll 1 \;, \end{array}\right. } \end{aligned}$$

(15)

which explains the concave and linear variation of $G_C$ as function of x/L observed in Fig. 4a and b, respectively. The growth of the bubbles along the tube is therefore a function of $n_P=\Delta P/P_L$.

Fig. 6

Variation of the bubble growth $G_\textrm{C}$, scaled with $n_P=\Delta P/P_L$, with x/L. Results are plotted for the same sets of $n_P$ for two different values of $P_L$. The left and right figures correspond to $\gamma =0$ and $0.03\,\mathrm{N/m}$, respectively. In each plot, the line corresponds to $P_L=1\,\textrm{kPa}$ and the symbols correspond to $P_L=100\,\textrm{kPa}$.

In Figure 6, we plot $G_\textrm{C}/n_P$ for the two outlet pressures $P_L$ with the same sets of values of $n_P$ for (a) $\gamma = 0$ and (b) $\gamma = 0.3\,\mathrm{N/m}$. The plots show that the results for the same values of $n_P$ follow the same curves, irrespective of the outlet pressures $P_L$. Furthermore, for the non-zero surface tension case in Fig. 6 (b), $G_\textrm{C}$ also shows a periodic oscillation along x when both the $n_P$ and $P_L$ are small, that is, for $P_L=1\,\textrm{kPa}$ and $n_P\le 1$. In addition, there is no data point for $n_P\le 0.3$ with $P_L=1\,\textrm{kPa}$, as the movement of the bubbles stopped due to high capillary barriers. This suggests the existence of an effective threshold pressure, below which there will be no flow through the tube. This threshold depends on both $\gamma$ and $P_L$, which we will explore more in the following section. We show the different characteristics of flow in the videos provided in electronic supplementary material.

3.3 Effective Rheology

Equations 1 and 2 resist analytical solutions even in the case when there is only a single compressible bubble in the tube. This is due to the pressure in the compressible bubble being inversely proportional to the difference in position of the two menisci surrounding it, whereas the motion of the two surrounding incompressible fluids is determined by the cosine of the positions of the same menisci. These equations, even in this simplest case, are therefore highly nonlinear with an essential singularity lurking in the very neighborhood where we seek solutions. We therefore stick to numerical analysis in the following.

Fig. 7

Plot of the flow rates for the total ($Q_\textrm{T}^\textrm{i,o}$), compressible ($Q_\textrm{C}^\textrm{i,o}$) and incompressible ($Q_\textrm{I}^\textrm{i,o}$) fluids at the inlet (left column) and at the outlet (right column) for $P_L=1\,\textrm{kPa}$ as a function of $\Delta P$. The different sets in each plot correspond to different values of the surface tension indicated in the legends. The quantities are divided with $Q^\prime =1\,\textrm{m}^3\mathrm{/s}$ and $P^\prime =1\,\textrm{kPa}$, respectively, to make them dimensionless. The dashed line in each plot has a slope 1.

Due to the volumetric growth of the compressible bubbles during their flow towards the outlet, the volumetric flow rate varies along the tube. In addition, this volumetric growth is a function of the pressures, making the average saturation and the effective viscosity of the two fluids inside the tube pressure dependent. These two mechanisms together control the effective rheological behavior of the steady-state flow. In Fig. 7, we show the variation of the volumetric flow rates ($Q_\textrm{T}^\textrm{i,o}$, $Q_\textrm{C}^\textrm{i,o}$, $Q_\textrm{I}^\textrm{i,o}$) as functions of the pressure drop $\Delta P$ for the outlet pressure $P_L=1\,\textrm{kPa}$ and for different values of the surface tension ($\gamma$). Note the differences between the inlet and outlet flow rates for the total and for the each component of flow. For the incompressible fluid, there is no increase in the outlet flow rate compared to its inlet flow rate (third row in Fig. 7) whereas there is a noticeable increase in the outlet flow rate of the compressible fluid (second row in Fig. 7). This increase in $Q_\textrm{C}^\textrm{o}$ effectively increases the total flow rate at the outlet (first row in Fig. 7). The dashed line in Fig. 7 has a slope equal to 1. The total flow rates show deviations from this dashed line. For the inlet, $Q_\textrm{T}^{i}$ shows small deviations from the dashed line for $\gamma >0$ at small $\Delta P$. Whereas at the outlet, the deviations are significantly higher due to the increase in the volumetric growth of the compressible fluid.

Another point to note in Fig. 7 is that there is a minimum value of $\Delta P$, below which there is no data point available. This is due to the existence of a threshold pressure below which the flow stops. In the supplementary material we show a simulation video in this regime where one can observe that the flow of the bubbles stops at a certain time step. The threshold is due to the capillary forces at the menisci between the two fluids that create capillary barriers at the narrowest points along the tube. Such threshold was also observed in the case of two-phase flow of two incompressible fluids in a tube with variable radius (Sinha et al. 2013). There, it was shown analytically that the average flow rate Q in the steady state varies with the applied pressure drop $\Delta P$ as, $Q\sim \sqrt{\Delta P^2 - P_\textrm{t}^2}$ where $P_\textrm{t}$ is the effective threshold pressure. When $|{\Delta P} | - P_\textrm{t} \ll P_\textrm{t}$, this relationship translates to $Q\sim \sqrt{|{\Delta P} | - P_\textrm{t}}$, that is, the flow rate varies with the excess pressure drop to the power of 0.5. The threshold pressure depends on the surface tension and on the configuration of the menisci positions inside the tube. If the total capillary barrier is higher than the applied pressure drop, the flow stops. This is similar here for the two-phase flow with one of the fluids being compressible.

Fig. 8

Plot of the volumetric inlet flow rate $Q_\textrm{T}^\textrm{i}$ as a function of the excess pressure drop $(\Delta P-P_\textrm{t})$ for $P_L=1\,\textrm{kPa}$ and $100\,\textrm{kPa}$, where the values of $P_\textrm{t}$ are obtained from a minimization of the least square fit error $\epsilon$. Here $Q^\prime =1\,\textrm{m}^3\mathrm{/s}$ and $P^\prime =1\,\textrm{kPa}$. The minimization is illustrated in the insets of a and b for $\gamma =0.05\,\mathrm{N/m}$ (green) and $0.09\,\mathrm{N/m}$ (purple). The dashed lines in a and b have a slope 1 whereas in c, the lower and upper dashed lines have slopes 1 and 1.3, respectively.

We assume a general relation between the average volumetric flow rates $Q_\textrm{T}^\textrm{i,o}$ and the pressure drop $\Delta P$ as,

$$\begin{aligned} \displaystyle Q_\textrm{T}^\textrm{i,o} \sim (\Delta P - P_\textrm{t})^{\beta _\textrm{i,o}} \end{aligned}$$

(16)

where $\beta _\textrm{i,o}$ is the corresponding exponent. In order to find both the effective threshold pressure $P_T$ and the exponent $\beta _\textrm{i,o}$ from the measurements of $Q_\textrm{T}^\textrm{i,o}$, we adopt an error minimization technique that was used in earlier studies (Sinha and Hansen 2012; Fyhn et al. 2021). There we choose a series of trial values for $P_\textrm{t}$ and calculate the mean square error $\epsilon$ for the linear least square fit by fitting the data points with $\log (Q)\sim \log (\Delta P-P_\textrm{t})$. Then we select the value of $P_\textrm{t}$ that corresponds to the minimum value of $\epsilon$, implying the best fit of the data points with Equation 16. This is illustrated in the insets of Fig. 8 (a) and (b). The slope for the selected threshold $P_\textrm{t}$ provides the exponent $\beta _\textrm{i,o}$. The variation of the total inlet and outlet flow rates $Q_\textrm{T}^\textrm{i,o}$ with the excess pressure drop $(\Delta P-P_\textrm{t})$ are plotted in Fig. 8 for the two outlet pressures $P_L=1$ and $100\,\textrm{kPa}$. The data sets show agreement with Equation 16 with the selected values of $P_\textrm{t}$ and $\beta$. There is a noticeable difference between the slopes for the inlet and outlet flow rates for $P_L=1\,\textrm{kPa}$ whereas for $P_L=100\,\textrm{kPa}$ they are similar. For $P_L=100\,\textrm{kPa}$ the data points for both $Q_\textrm{T}^\textrm{i}$ and $Q_\textrm{T}^\textrm{o}$ follow a slope of $\approx 1.0$ whereas for $P_L=1\,\textrm{kPa}$, the data points for $Q_\textrm{T}^\textrm{i}$ and $Q_\textrm{T}^\textrm{o}$ follow the slopes of $\approx 1.0$ and 1.3, respectively. These are indicated by the dashed lines in the figures.

Fig. 9

Variation of the threshold pressure $P_\textrm{t}$ and the exponents $\beta _\textrm{i,o}$ as functions of the effective surface tension $\gamma$ for $P_L=1$ and $100\,\textrm{kPa}$. $P_\textrm{t}$ increases with the increase of $\gamma$ and the values are much higher for $P_L=1\,\textrm{kPa}$ compared to $P_L=100\,\textrm{kPa}$. The exponent $\beta _\textrm{i}$ for the inlet flow rate are close to 1 for both the values of $P_L$ whereas for the outlet flow rate $\beta _\textrm{o}\approx 1.3$ for $P_L=1\,\textrm{kPa}$. For $P_L=100\,\textrm{kPa}$, $\beta _\textrm{o}$ remains close to $\beta _\textrm{i}$. The dashed horizontal lines indicate the value 1.0 of the y axis.

The variations of $P_\textrm{t}$ and $\beta _\textrm{i,o}$ with the surface tension $\gamma$ are plotted in Fig. 9. The data points were calculated by considering different ranges of $\Delta P$ and taking averages over the ranges, and the corresponding standard deviations are plotted as error bars. The threshold pressure $P_\textrm{t}$ is zero at $\gamma =0$ and then increases gradually with $\gamma$ which shows that the threshold appears due to capillary forces. The increase in $P_\textrm{t}$ with $\gamma$ appears to be linear here which is similar to the case of two incompressible fluids, where the linear dependence of $P_\textrm{t}$ on the surface tension was shown analytically (Sinha et al. 2013). Additionally for the compressible flow here, the thresholds also depend on the outlet pressure $P_L$. For the lower outlet pressure $P_L=1\,\textrm{kPa}$, the thresholds are systematically higher compared to those for $P_L=100\,\textrm{kPa}$ for the whole range of $\gamma$. Furthermore, the exponents $\beta _\textrm{i,o}$ also depend on the outlet pressure as seen from Figs. 9 (b) and (c). The difference is more visible for the exponents related to the outlet flow rates than the inlet. For the inlet flow rate, $\beta _\textrm{i}$ has values around $\approx 0.95$ and 1.02 for $P_L=1\,\textrm{kPa}$ and $100\,\textrm{kPa}$, respectively, showing almost linear dependence for both the cases. For the outlet flow rates, $\beta _\textrm{o}$ remains close to $\beta _\textrm{i}$ for $P_L=100\,\textrm{kPa}$ whereas for $P_L=1\,\textrm{kPa}$, $\beta _\textrm{o}$ increases to $\approx 1.3$. This increase in $\beta _\textrm{o}$ compared to $\beta _\textrm{i}$ reflects the dependence of the volumetric growth $G_\textrm{C}(x)$ of the bubbles on $P_L$, indicating an underlying dependence of the rheological behavior on the absolute inlet or outlet pressures. However, at this point we are unable to describe how the two parameters $P_\textrm{t}$ and $\beta$ scale with $P_L$, which needs further study. In addition, we have only considered an intermediate volumetric fractional flow $F_\textrm{C} = 0.4$ here, which also controls $\beta$ and $P_\textrm{t}$ (Roy et al. 2021). If the fractional flow or the saturation is made near to either 0 or 1, the system will approach single-phase flow and the linearity in the rheology of the Newtonian fluids should be retrieved.

Compared to the study of a single capillary tube here, a porous medium is composed of many interconnected pores of different sizes. Existing studies of two-phase flow in porous media have shown the existence of different power-law regimes for the relation between volumetric flow rate and pressure drop. These regimes are characterized by different exponents. The studies involve experiments (Tallakstad et al. 2009, 2009; Rassi et al. 2011; Sinha et al. 2017; Gao et al. 2020; Zhang et al. 2021), Lattice Boltzmann simulations (Yiotis et al. 2013), pore-network modeling (Sinha and Hansen 2012; Sinha et al. 2017) and analytical calculations (Tallakstad et al. 2009; Sinha and Hansen 2012; Roy et al. 2019; Zhang et al. 2021). There are three regimes, an intermediate nonlinear regime where the flow rate Q increases at a rate much faster than the applied pressure drop $\Delta P$ with a power law exponent larger than one and up to around 2.5. There are in addition two linear regimes for either smaller (Yiotis et al. 2013; Gao et al. 2020; Zhang et al. 2021) or larger (Yiotis et al. 2013; Sinha and Hansen 2012; Sinha et al. 2017) volumetric flow rates than the nonlinear regime. This allows the definition of a lower and upper crossover pressure drop. The origin of the power law in a porous network and the crossovers to different regimes, can be explained by two dominant factors, the rheology of individual pores and the distribution of the threshold pressures in the network (Roy et al. 2019). A simple explanation can be drawn from a disordered network of threshold resistors (Roux and Herrmann 1987) where each resistor has a threshold voltage to start conducting the current. In a network with links with a distribution of thresholds, there will be a regime when new conducting paths will appear while increasing the global pressure drop. The increase in the flow rate through each path together with the increase in the number of paths leads to an effective increase in Q faster than $\Delta P$. This results in the nonlinear exponent being higher than 1, the value of which depends on the distribution of the thresholds in each link (Roy et al. 2019). The linear regime above this nonlinear regime appears from all the available paths being conducting whereas the linear regime below appears from the flow being in single percolating channels, which are governed by the rheology of individual pores. According to this explanation, the experimental (Gao et al. 2020; Zhang et al. 2021) and numerical (Yiotis et al. 2013) observations of two-phase flow in porous media showing linear variation of flow rate in the low pressure regime therefore indicate that the flow in the single channels consisting of many pores are linear, which is similar to what we have found for the lower outlet pressure in the present compressible/incompressible flow case.

4 Conclusions

We have studied the flow of alternating compressible bubbles and incompressible droplets through a capillary tube with variable radius. The motion of the bubbles was given by the model Equations 1 and 2, thus assuming the compressible fluid to be an ideal gas with zero viscosity, whereas the incompressible fluid is Newtonian. The incompressible fluid is more wetting than the compressible gas, but not to a degree that films form. We switch between injecting the compressible and incompressible fluid at intervals so that the fractional flow rate is essentially constant at the inlet. We fix the pressure drop along the tube in addition to an ambient pressure. This creates steady-state flow conditions in the tube.

The compressible bubbles expand as they move from the higher pressure region at the inlet towards the lower pressure at the outlet. This expansion accelerates the incompressible fluid, thus making the volumetric flow rate larger at the outlet than at the inlet. The lower the ambient pressure is, the stronger this effect is.

We measure volumetric flow rate at the inlet, finding essentially a linear relationship between the volumetric flow rate and the pressure drop. However, there is a threshold pressure that needs to be overcome in order to have flow through the tube. At the outlet, we find that the volumetric flow rate is still linear in the excess pressure drop when the ambient pressure is low. However, when the ambient pressure is high, the volumetric flow rate at the outlet becomes proportional to the excess pressure to a power of around 1.3.

This behavior is very different from that of two incompressible fluids moving through a corresponding tube: Here the volumetric flow rate, being the same at the inlet and the outlet, is proportional to the square root of the excess pressure.

We expected the flow rate-pressure drop constitutive relations to be different in this compressible/incompressible case than that of two incompressible fluids. However, that we should find linearity was a big surprise. A precise explanation as to why this is so, is still lacking.

Besides these surprising results, this work makes a first step in modeling of compressible/incompressible fluid mixtures in dynamic network models. We may then envision using more sophisticated equations of state for the compressible fluid beyond the ideal gas law. This allows the consideration of, e.g., phase transitions such as boiling and condensation in porous media.

Finally, we like to point out that we have not considered film-flow or contact line pinning in this study. Pinning of the contact lines will change the relation between the capillary force at the interfaces and the shape of the tube, which is determined by two more fixed parameters: the surface tension and the wetting angle. If the pinning is due to surface roughness at scales much finer than the variations in the tube radius, this will result in an effective wetting angle different from the one expected for smooth surfaces (see Blunt 2017, pp 11-14). Hence, we do not expect our results to be qualitatively different in this case. If, on the other hand, the roughness is on the same scale as the radius variations, we are dealing with a tube that essentially has a different shape than the one we are considering. Even in this case, we do not expect qualitative changes from the results we report. However, there will most probability be quantitative changes. Different tube shapes were studied by Lanza et al. for an immiscible mixture of a yield stress fluid and a Newtonian fluid, finding a quantitative difference between different tube shapes (Lanza et al. 2022; Talon et al. 2014). However, both the fluids were incompressible in that study, therefore a possible future extension of both of these studies would be to include a compressible fluid together with a non-Newtonian fluid in a capillary.

The film flow on the other hand will introduce parallel components of the two fluids in the system whereas in our present problem the two fluid components are always in series combination. Depending on the thickness of the films, this may change the effective relationship between the flow rate and the pressure drop. Film flow can be observed when the pores contain rough grain surfaces and corners (Chen et al. 2018; Cejas et al. 2018) or when a fluid phase is completely wetting (Aursjø et al. 2014), whereas in case of drainage dominated flow the film flow may be neglected. Experiments have also shown that gravity plays a role in controlling the active zone of film flow in a porous media (Moura et al. 2019). Experiments with the same porous media with different types of fluids have shown that the nonlinear exponent $\beta$ was smaller for the fluids that show strong film flow (Aursjø et al. 2014) compared to those without film flow (Tallakstad et al. 2009). Fluid wettabilities in this context strongly affects the appearance of films as well as the rheological nonlinearity in general (Zhang et al. 2022; Fyhn et al. 2021). How the introduction of films or changing the wettability of the fluids will affect the results of the present study is therefore a question for the future.