3.1. Geometry of a trapped cluster

When describing a cluster’s geometry, it is convenient to label the obstacles along its perimeter in counterclockwise order as shown in figure 2, such that the $n$ th opening of the cluster lies between the $n$ th and the $(n+1)$ th obstacle. We equip each opening with a local coordinate system $\{ \boldsymbol{\hat {t}}_n, \boldsymbol{\hat {n}}_n \}$ , where $\boldsymbol{\hat {t}}_n$ is the unit vector pointing from the centre of the $n$ th obstacle, $\boldsymbol{r}_{\textit{obs}, n}$ , towards the centre of the $(n+1)$ th obstacle, $\boldsymbol{r}_{\textit{obs},n+1}$ , and $\boldsymbol{\hat {n}}_n$ is found by rotating $\boldsymbol{\hat {t}}_n$ $90^\circ$ counterclockwise. The orientation of the $n$ th local coordinate system relative to the global $\{ \boldsymbol{\hat {x}}, \boldsymbol{\hat {y}}\}$ coordinate system is encoded in the orientation angle $\xi _n$ between the global $x$ -axis and $\boldsymbol{\hat {t}}_n$ . We emphasise that all these quantities are independent of ${\textit{Bo}}$ and can be calculated once and for all before solving (3.1).

Now, assuming the fluid interfaces to be circular, each of them can be described by their radius, $R_n$ , and their centre, $\boldsymbol{r}_{c,n}$ – see figure 3. Since a circular arc has constant curvature, we note that this assumption neglects the possibility that the pressure jump across the interface may in fact vary along the interface. This effect is expected to be largest for a cluster opening which is parallel to the mean flow direction, and for such an opening, one can show that the difference between the curvature at the bottom and the top of the interface is roughly proportional to the Bond number, the relative difference of the fluids’ densities as well as the width of the opening. As a consequence, the assumption of circular interfaces is only valid as long as the product of these quantities is small. One way of calculating the quantities $R_n$ and $\boldsymbol{r}_{c,n}$ in terms of the opening angle $\omega _n$ is to parametrise the curves along the circumference of the $n$ th obstacle and the $n$ th interface; the two curves should intersect when the angles of the parametrisations are $\omega _n$ and $\varOmega _n$ , respectively, and at the intersection point, the angle between the two curves should be the contact angle, $\theta _0$ . In the local coordinate system, a parametrisation of the circumference of the $n$ th obstacle is given by

(3.2) \begin{align} \boldsymbol{r}(\omega ) = \boldsymbol{r}_{\textit{obs},n} + r \cos (\omega ) \boldsymbol{\hat {t}}_n + r \sin (\omega ) \boldsymbol{\hat {n}}_n, \end{align}

while the $n$ th fluid interface can be parametrised as

(3.3) \begin{align} \boldsymbol{r}(\varOmega ) = \boldsymbol{r}_{c,n} - R_n \cos (\varOmega ) \boldsymbol{\hat {t}}_n + R_n \sin (\varOmega ) \boldsymbol{\hat {n}}_n. \end{align}

We note that the obstacle parametrisation runs counterclockwise and that the interface parametrisation runs clockwise, in accordance with the signs of $\omega _n$ and $\varOmega _n$ defined in figure 3. The angle $\varOmega _n$ is found by requiring that the angle between the two curves’ tangent vectors $\boldsymbol{\hat {t}}_{\textit{obs},n}$ and $\boldsymbol{\hat {t}}_{\textit{int},n}$ should equal the contact angle, and a straightforward calculation shows that

(3.4) \begin{align} \varOmega _n = \theta _0 - \omega _n. \end{align}

To calculate the radius of the fluid interface, we make use of the fact that all obstacles have the same radius. This implies that the centre of the interface must lie half-way between the obstacles, that is, on the perpendicular bisector of line connecting the two obstacle centres, and hence that $(\boldsymbol{r}_{c,n} - \boldsymbol{r}_{\textit{obs},n})\boldsymbol{\cdot }\boldsymbol{\hat {t}}_n = r + d_n/2$ . Equating (3.2) and (3.3) at the point where $\omega = \omega _n$ and $\varOmega = \varOmega _n$ , and taking the dot product with $\boldsymbol{\hat {t}}_n$ therefore implies that the normalised interface radius is given by the expression

(3.5) \begin{align} \frac {R_n}{r} = \frac {d_n/r + 2\big ( 1 - \cos (\omega _n) \big ) }{2 \cos (\varOmega _n)}. \end{align}

The geometric description of the interface is completed by equating (3.2) and (3.3) at the intersection point and taking the dot product with $\boldsymbol{\hat {n}}_n$ ; when combined with the above-mentioned symmetry argument, this operation yields

(3.6) \begin{align} \frac {\boldsymbol{r}_{c,n}}{r} = \frac {\boldsymbol{r}_{\textit{obs},n}}{r} + \bigg (1 + \frac {d_n}{2 r} \bigg ) \boldsymbol{\hat {t}}_n - \bigg ( \frac {R_n}{r}\sin (\varOmega _n) - \sin (\omega _n) \bigg ) \boldsymbol{\hat {n}}_n. \end{align}

3.2. Area of a trapped cluster

To derive the area of a trapped cluster, we again consider the cluster depicted in figure 2. Using Green’s theorem, its area, $A_{\textit{cluster}}$ , can be calculated as the line integral

(3.7) \begin{align} A_{\textit{cluster}} = \frac {1}{2} \oint _{\mathcal{C}} \big ( -y \, \mathrm{d} x + x \, \mathrm{d} y\big ), \end{align}

where $\mathcal{C}$ is the closed curve along the cluster’s circumference oriented counterclockwise. We note that the integral only equals the cluster area in the case where the cluster contains no embedded obstacles, that is, obstacles that are completely surrounded by wetting fluid, which may be present in the case of large clusters. However, since our only application of the cluster area is to make sure that it remains constant when the Bond number changes, the constant contributions to the area from embedded obstacles can be neglected for the present purpose. To make progress from here, we decompose $\mathcal{C}$ into contributions from the fluid interfaces and the obstacles, as shown in figure 2, which gives that

(3.8) \begin{align} A_{\textit{cluster}} = \frac {1}{2} \sum _{n = 1}^{N} \bigg ( \oint _{\mathcal{C}_{\textit{int},n}} \big ( -y \, \mathrm{d} x + x \, \mathrm{d} y \big ) + \oint _{\mathcal{C}_{\textit{obs},n}} \big ( -y \, \mathrm{d} x + x \, \mathrm{d} y \big ) \bigg ). \end{align}

To calculate the contribution from the $n$ th fluid interface, we employ the parametrisation (3.3) with $\varOmega _n \leqslant \varOmega \leqslant \pi - \varOmega _n$ . Carrying out the integral (see Appendix A for a detailed derivation) yields

(3.9) \begin{align} \int _{\mathcal{C}_{\textit{int},n}} \big ( -y \, \mathrm{d} x + x \, \mathrm{d} y \big ) = - 2 R_n \cos \big (\varOmega _n \big ) \big ( \boldsymbol{r}_{\text{c},n} \boldsymbol{\cdot }\boldsymbol{\hat {n}}_n \big ) - R_n^2 \big ( \pi - 2\varOmega _n \big ). \end{align}

For the integral along the obstacle circumference, instead of using the local coordinate system, it turns out to be beneficial to use the global $\{ \boldsymbol{\hat {x}}, \boldsymbol{\hat {y}}\}$ -coordinate system to parametrise the curve segment. One possible choice is

(3.10) \begin{align} \boldsymbol{r}(s) = \boldsymbol{r}_{\textit{obs},n} - r \cos (s) \boldsymbol{\hat {x}} + r \sin (s) \boldsymbol{\hat {y}}, \quad \text{where} \quad s_{1, n} \leqslant s \leqslant s_{2, n}, \end{align}

in which the upper limit of the parametrisation is $s_{2,n} = \pi - \xi _n - \omega _n$ , while the lower limit is $s_{1,n} = 2\pi - \xi _{n-1} + \omega _{n-1}$ provided that this value is smaller than $s_{2,n}$ ; if it is larger than $s_{2,n}$ , $2\pi$ should be subtracted from $s_{1,n}$ . This technicality only matters when evaluating the term $r^2 (s_{2,n} - s_{1,n})$ in (3.11), and we note that the subtraction of $2\pi$ is in fact irrelevant when establishing (3.1), since it disappears when differentiating with respect to ${\textit{Bo}}$ . Inserting the parametrisation (3.10) into the line integral, the contribution from the $n$ th obstacle segment can be calculated in exactly the same way as the integral (3.9) and turns out to be

(3.11) \begin{align} \begin{split} \int _{\mathcal{C}_{\textit{obs}},n} \big ( -y \, \mathrm{d} x + x \, \mathrm{d} y \big ) & = r y_{\textit{obs},n} \big ( \cos (s_{2,n}) - \cos (s_{1,n}) \big ) \\ & + r x_{\textit{obs},n} \big ( \sin (s_{2,n}) - \sin (s_{1,n}) \big ) \\ & - r^2 (s_{2,n} - s_{1,n}), \end{split} \end{align}

where $x_{\textit{obs},n}$ and $y_{\textit{obs},n}$ are the $x$ - and $y$ -coordinates of $\boldsymbol{r}_{\textit{obs},n}$ , respectively. In combination, the results (3.8), (3.9) and (3.11) imply that the cluster area must satisfy the relation

(3.12) \begin{align} \begin{split} 2A_{\textit{cluster}}/r^2 = \sum _{n = 1}^{N} f_n (\omega _n, \omega _{n-1}), \end{split} \end{align}

in which the function $f_n(\omega _n, \omega _{n-1})$ is defined as

(3.13) \begin{align} \begin{split} f_n(\omega _n, \omega _{n-1}) \equiv & - \frac {2 R_n}{r} \cos \big (\varOmega _n \big ) \big ( (\boldsymbol{r}_{c, n} / r) \boldsymbol{\cdot }\boldsymbol{n}_n \big ) - \bigg (\frac {R_n}{r} \bigg )^2 \big ( \pi - 2\varOmega _n \big ) \\ & + \frac {y_{\textit{obs},n}}{r} \big ( \cos (s_{2,n}) - \cos (s_{1,n}) \big ) \\ & + \frac {x_{\textit{obs},n}}{r} \big ( \sin (s_{2,n}) - \sin (s_{1,n}) \big ) \\ & - (s_{2,n} - s_{1,n}). \end{split} \end{align}

This completes the calculation of the cluster area.

3.3. Pressure difference across a trapped cluster

As long as a cluster is trapped, the pressure just outside its $n$ th opening is related to the pressure just outside its $i$ th opening. To show this, we first define the outer and inner points of the $n$ th fluid interface as $\boldsymbol{r}_{\textit{out},n} \equiv \boldsymbol{r}_{\textit{int},n} - \delta \boldsymbol{\hat {n}}_n$ and $\boldsymbol{r}_{\textit{in},n} \equiv \boldsymbol{r}_{\textit{int},n} + \delta \boldsymbol{\hat {n}}_n$ , and equivalently for the $i$ th interface. Here, $\boldsymbol{r}_{\textit{int},n} \equiv \boldsymbol{r}_{c,n} + R_n \boldsymbol{\hat {n}}_n$ is the midpoint of the $n$ th interface as illustrated in figure 3, and $\delta$ is a length scale which is very large relative to the width of the interface and very small relative to the obstacle radius. Next, we consider the momentum balance across the $n$ th interface, which requires that the equation

(3.14) \begin{align} \begin{split} \big ( p_{\textit{out},n} - p_{{in}, n} \big ) \boldsymbol{\hat {n}}_n = & \gamma \kappa _n \boldsymbol{\hat {n}}_n + \mu _{{nw}} \Big (\boldsymbol{\nabla }\boldsymbol{u}_{\textit{out},n} + (\boldsymbol{\nabla }\boldsymbol{u}_{\textit{out},n})^\top \Big ) \boldsymbol{\hat {n}}_n \\ & - \mu _{{w}} \Big (\boldsymbol{\nabla }\boldsymbol{u}_{\textit{in},n} + (\boldsymbol{\nabla }\boldsymbol{u}_{\textit{in},n})^\top \Big ) \boldsymbol{\hat {n}}_n \end{split} \end{align}

is satisfied (see e.g. the derivation in § 1.1.2 of the book by Gross & Reusken Reference Gross and Reusken2011). In this equation, $\kappa _n \equiv 1/R_n$ is the curvature of the $n$ th interface, $\boldsymbol{u}$ is the velocity field, and the subscripts ‘out, $n$ ’ and ‘in, $n$ ’ signify that the sub-scripted quantity is evaluated at $\boldsymbol{r}_{{out}, n}$ and $\boldsymbol{r}_{{in}, n}$ , respectively. Since $\boldsymbol{u}$ is continuous across the interface, we note that the fluid inside the trapped cluster may in fact circulate. From the simulations presented in §§ 4 and 5, we have, however, found that the fluid velocity and its gradient inside the trapped cluster are always very small relative to the corresponding quantities of the outside flow, and for that reason, we neglect the last term on the right-hand side of (3.14). Interestingly, this implies that the viscosity of the wetting fluid does not impact the stability of the cluster and, as a consequence, that the ratio of the two fluids’ viscosities cannot be a relevant parameter for cluster stability (provided, as in this case, that substantial circulation inside the cluster can be ruled out). In addition, a crude order-of-magnitude estimate shows that $\mu _{{nw}} |\boldsymbol{\nabla }\boldsymbol{u}_{\textit{out},n} | /(\gamma \kappa _n) = O({\textit{Ca}})$ , where ${\textit{Ca}} \equiv \mu _{{nw}} v_0 / \gamma$ is the capillary number of the flow in the non-wetting fluid and $v_0$ is a typical velocity of the flow. Hence, as long as the capillary number is small, the second term on the right-hand side of (3.14) can be neglected and the pressure differences across the $n$ th and $i$ th interfaces are given by the well-known Young–Laplace equation, i.e.

(3.15) \begin{align} \left \{ \begin{aligned} p_{\textit{out},n} - p_{\textit{in},n} & = \kappa _n \gamma , \\[3pt] p_{\textit{out},i} - p_{\textit{in},i} & = \kappa _i \gamma . \end{aligned} \right . \end{align}

As is evident from these equations, we choose the curvature to be positive when the pressure on the outside is larger than on the inside, that is, $\kappa _n$ is positive whenever the non-wetting side of the interface is convex. The final step in relating $p_{\textit{out},n}$ to $p_{\textit{out},i}$ consists in noting that the pressure inside the trapped cluster must be hydrostatic; the contribution from viscous effects to the pressure difference is insignificant because the fluid – to a very good approximation – lies still there and, hence,

(3.16) \begin{align} p_{\textit{in},n} - p_{\textit{in},i} = \rho _{{w}} g (y_{\textit{int}, n} - y_{\textit{int,i}}), \end{align}

and if we subtract the two equations in (3.15), use (3.16) to eliminate the pressure difference inside the cluster and normalise by $\rho _{{w}} g r$ , we find that

(3.17) \begin{align} \frac {p_{\textit{out},n} - p_{\textit{out},i}}{\rho _{{w}} g r} = \frac {\kappa _n r - \kappa _i r}{{\textit{Bo}}} + \frac {y_{\textit{int},n} - y_{\textit{int},i}}{r}. \end{align}

At a first glance, one would think that this relation imposes $N (N - 1)/2$ constraints on the geometry of the cluster, because it should hold for all possible combinations of $n$ and $i$ . However, if it holds for $i = 1$ and all $n \geqslant 2$ , then the relation is automatically satisfied for all possible choices of $n$ and $i$ . This is implied by the linearity of the relation. As such, (3.17) only imposes $N-1$ constraints on the cluster geometry.

In what follows, we will assume that the left-hand side of (3.17) is independent of the Bond number. As we envision changing the Bond number by varying either $g$ or $\gamma$ , this is reasonable when recalling that a creeping flow is accurately described by the Stokes equation and therefore linear in the driving force and otherwise determined by the geometry. Since the pressure difference is normalised by $g$ , changing the value of this parameter does not have a direct effect on the left-hand side. Moreover, as we only consider the range ${\textit{Bo}} \leqslant {\textit{Bo}}_{{crit}}$ , the changes in the cluster geometry are necessarily small. In particular, no pathways for the non-wetting fluid are created nor destroyed in this range. As such, the geometrical contribution to the change in the normalised pressure difference is expected to be small, which it indeed is. Figure 4 shows an example of this by comparing a trapped clusters in two different states in which it is either very far or very close to breaking. The validation of the quasistatic model against numerical simulations in § 4 provides further evidence that the assumption is indeed satisfied to a good approximation.

Figure 4.

Normalised pressure field around a trapped cluster for ${\textit{Bo}} = 0.025$ and ${\textit{Bo}} = 0.125$ , corresponding to capillary numbers ${\textit{Ca}} = 0.00025$ and ${\textit{Ca}} = 0.00125$ , respectively. In the two cases, the cluster is either very far (panel a) of very close (panel b) to breaking, yet the pressure fields around the cluster are essentially the same. The pressure inside the cluster is larger in the small- ${\textit{Bo}}$ case than in the large- ${\textit{Bo}}$ case, because the Bond number is varied by varying the surface tension, while keeping the gravitational acceleration constant. The pressure field is calculated using the DNS model described in § 4.1 using a circular cluster with radius $r_{{circle}}/r = 4$ as the initial condition.

3.4. System of differential equations governing the cluster geometry

Equations (3.12) and (3.17) govern the geometry of the cluster as the Bond number changes. As we have shown, all geometric changes are ultimately determined by changes in the opening angles, $\boldsymbol{\omega } = [\omega _1, \omega _2, \ldots , \omega _N]$ , which in turn depend on a long list of parameters. However, as we are presently only concerned with changes in the opening angles induced by changes in ${\textit{Bo}}$ , we consider $\omega _n$ a function of this variable only. Hence, if we require the mass of the cluster to be independent of the Bond number, taking the derivative of (3.12) with respect to ${\textit{Bo}}$ , yields

(3.18) \begin{align} 0 = \sum _{n = 1}^{N} \bigg ( \frac {\partial f_n}{\partial \omega _n} + \frac {\partial f_{n+1}}{\partial \omega _n} \bigg ) \frac {\partial \omega _n}{\partial ({\textit{Bo}})}. \end{align}

Here, periodic indexing is employed such that $f_{N+1}(\omega _{N+1}, \omega _{N}) = f_1(\omega _1, \omega _{N})$ . This is the first of the $N$ equations that constitute (3.1); the remaining $N-1$ equations are derived from the relation (3.17). Setting $i = 1$ and making the assumption that the left-hand side is independent of the Bond number, taking the derivative with respect to ${\textit{Bo}}$ yields

(3.19) \begin{align} \begin{split} 0 & = \frac {\mathrm{d}}{\mathrm{d} ({\textit{Bo}})} \bigg ( \frac {p_{\textit{out},n} - p_{\textit{out},1}}{\rho _{{w}} g r} \bigg ) \\[5pt] & = \frac {\mathrm{d}}{\mathrm{d} ({\textit{Bo}})} \bigg ( \frac {\kappa _n r - \kappa _1 r}{{\textit{Bo}}} + \frac {y_{\textit{int},n} - y_{\textit{int},1}}{r} \bigg ) \\[5pt] & = \frac {1}{{\textit{Bo}}} \frac {\mathrm{d} (\kappa _n r)}{\mathrm{d} {\textit{Bo}}} - \frac {1}{{\textit{Bo}}} \frac {\mathrm{d} (\kappa _1 r)}{\mathrm{d} {\textit{Bo}}} - \frac {\kappa _n r - \kappa _1 r}{{\textit{Bo}}^2} + \frac {\mathrm{d} (y_{\textit{int},n}/r)}{\mathrm{d} ({\textit{Bo}})} - \frac {\mathrm{d} (y_{\textit{int},1}/r)}{\mathrm{d} ({\textit{Bo}})} \end{split} \end{align}

for $n = 2, 3, \ldots , N$ . Rearranging terms then leads to the equation

(3.20) \begin{align} \bigg ( \frac {\mathrm{d} (y_{\textit{int},1}/r)}{\mathrm{d} \omega _1} + \frac {1}{{\textit{Bo}}} \frac {\mathrm{d} (\kappa _1 r)}{\mathrm{d} \omega _1} \bigg ) \frac {\partial \omega _1}{\partial ({\textit{Bo}})} - \bigg ( \frac {\mathrm{d} (y_{\textit{int},n}/r)}{\mathrm{d} \omega _n} + \frac {1}{{\textit{Bo}}} \frac {\mathrm{d} (\kappa _n r)}{\mathrm{d} \omega _n} \bigg ) \frac {\partial \omega _n}{\partial ({\textit{Bo}})} = \frac {\kappa _1 r - \kappa _n r}{({\textit{Bo}})^2}, \end{align}

which, in combination with (3.18), constitutes the system of $N$ coupled equations (3.1) governing the evolution of the opening angles as a function of the Bond number. To arrive at the matrix formulation (3.1), we simply define

(3.21a) \begin{align} a_n & \equiv \frac {\partial f_n}{\partial \omega _n} + \frac {\partial f_{n+1}}{\partial \omega _n}, \\[-12pt]\nonumber \end{align}

(3.21b) \begin{align} b_n & \equiv \frac {\mathrm{d} (y_{\textit{int},n}/r)}{\mathrm{d} \omega _n} + \frac {1}{{\textit{Bo}}}\frac {\mathrm{d} (\kappa _n r)}{\mathrm{d} \omega _n}, \\[-12pt]\nonumber \end{align}

(3.21c) \begin{align} c_n & \equiv \frac {\kappa _1 r - \kappa _n r }{({\textit{Bo}})^2}, \end{align}

and set the matrix $\unicode{x1D63C}$ and vector $\boldsymbol{c}$ to

(3.22) \begin{align}\unicode{x1D63C} = \begin{bmatrix} a_1 & a_2 & a_3 & {\cdots} & a_{N-1} & a_{N} \\[3pt] b_1 & -b_2 & 0 & {\cdots} & 0 & 0 \\[3pt] b_1 & 0 & -b_3 & {\cdots} & 0 & 0 \\[3pt] \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\[3pt] b_1 & 0 & 0 & {\cdots} & -b_{N-1} & 0 \\[3pt] b_1 & 0 & 0 & {\cdots} & 0 & -b_{N} \\ \end{bmatrix} \!, \quad \boldsymbol{c} = \begin{bmatrix} 0 \\[3pt] c_2 \\[3pt] c_3 \\[3pt] \vdots \\[3pt] c_{N-1} \\[3pt] c_{N} \end{bmatrix}\!. \end{align}

Given an initial configuration, $(\boldsymbol{\omega }_0, {\textit{Bo}}_0)$ , the set of equations (3.1) is easily integrated numerically with respect to the Bond number by any standard ordinary differential equation solver. As such, it only remains for us to clarify up to which Bond number one should integrate to predict the point of breaking. This is done in the following.

3.5. When does a cluster become unstable?

As long as the Bond number is smaller than its critical value, the interfaces of the cluster can balance out the pressure difference across them by adjusting their shapes. More precisely, if the pressure difference across the $n$ th interface is $\Delta p_n$ , then the interface will adjust its curvature such that the Young–Laplace equation is satisfied. Using the expression (3.5) for the radius of the interface, the expression (3.4) for the angle $\varOmega _n$ and normalising by $\rho _{{w}} g r$ , the Young–Laplace equation can be written as

(3.23) \begin{align} \frac {\Delta p_n}{\rho _{{w}} g r} = \frac {1}{{\textit{Bo}}} \, \frac {2 \cos (\theta _0 - \omega _n)}{d_n/r + 2\big ( 1 - \cos (\omega _n) \big )}. \end{align}

This is the (normalised) pressure difference that the $n$ th interface can sustain by opening up to the angle $\omega _n$ . Plots of this pressure difference as a function of $\omega _n$ are shown in figure 5 for $\theta _0 = \pi /3$ for different values of $d_n/r$ , and from this, it is clear that the pressure difference has a global maximum. One may interpret this result in the following way: if the cluster is initially stable, it can respond to an increasing pressure on its outer side, a decreasing pressure on its inner side or a weakening of the surface tension by moving inwards. Since the contact angle is constant, the new position further inwards will ensure that the curvature of the interface is larger and, therefore, that the interface becomes stronger. This, however, only works up to a certain critical angle, which we denote by $\omega _{\textit{crit},n}$ . Beyond this angle, further inwards movement will result in a smaller curvature and hence a weakening of the interface. Differentiating (3.23) with respect to $\omega _n$ , it is straightforward to show that the critical angle is given by the expression

(3.24) \begin{align} \omega _{\textit{crit},n} = \theta _0 - \sin ^{-1} \bigg ( \frac {2 \sin (\theta _0)}{d_n/r + 2} \bigg ). \end{align}

When the flow conditions are such that one of the cluster’s openings is pushed to its critical value, the cluster becomes unstable. If some small perturbation pushes this interface just a tiny bit further in, it will become weaker and therefore keep yielding to the outside pressure until it finally collapses. The critical value of the Bond number for a given cluster can thus be computed by solving (3.1) up to the point where $\omega _n = \omega _{\textit{crit},n}$ for some $n$ .

We note that this instability is the well-known Haines jump, discovered by Haines (Reference Haines1930) now nearly a century ago. Cieplak & Robbins (Reference Cieplak and Robbins1988) refer to this event as a ‘burst’ and, as they have shown, there are other ways in which an interface can be unstable. For example, an interface may split into two new interfaces if it touches another obstacle before reaching its critical angle. The detection of such an event could, in principle, be built into the numerical solution of the model proposed here. We have, however, not done so because the continuation after the touch would most likely require a full (highly expensive) recomputation of the outside pressure field. While the small movements of an interface do not lead to large changes in the surrounding pressure, a reconfiguration of the interfaces oftentimes does. For small clusters, we note that this simplification does not affect the computation of ${\textit{Bo}}_{{crit}}$ , because these do not have any embedded obstacles such that a touching event could take place. For large clusters, however, we acknowledge that the simplification may add to the uncertainty of the computed values of ${\textit{Bo}}_{{crit}}$ .

Figure 5.

Pressure difference across the $n$ th cluster opening given by (3.23) as a function of the angle $\omega _n$ for $\theta = \pi /3$ and different values of the non-dimensional inter-obstacle distance $d_n/r$ .

3.6. Varying other parameters than the Bond number

The Bond number is not the only parameter on which the trapped cluster’s geometry depends; for example, it also depends on the ratio of fluid densities, $\rho _{{w}}/\rho _{{nw}}$ , and the contact angle, $\theta _0$ , and using an approach similar to that mentioned earlier, it is possible to derive equations corresponding to (3.1) for predicting the change in $\boldsymbol{\omega }$ when these parameters are varied. In Appendix B, we show and validate that when the ratio $\rho _{{w}}/\rho _{{nw}}$ is varied while keeping all other variables fixed, the opening angles satisfy the system of equations

(3.25) \begin{align}\unicode{x1D63C} \frac {\partial \boldsymbol{\omega }}{\partial \left(\frac {\rho _{{w}}}{\rho _{{nw}}}\right)} = \boldsymbol{c}_{\rho _{{w}}/\rho _{{nw}}}, \end{align}

in which $\unicode{x1D63C}$ is the same as in (3.1) and $\boldsymbol{c}_{\rho _{{w}}/\rho _{{nw}}}$ is given by (B3). We note that this set of equations allows for cost-efficient computation of the critical Bond number as a function of the density ratio. Given a cluster geometry for particular values of ${\textit{Bo}}$ and $\rho _{{w}}/\rho _{{nw}}$ , one can first integrate $\boldsymbol{\omega }$ with respect to the density ratio up to the desired value using (3.25) before computing ${\textit{Bo}}_{{crit}}$ using (3.1). A set of equations that allows for changes in $\theta _0$ can likewise be derived, but we have not done so in the present work. A quasistatic model for the change in cluster geometry with respect to the parameters of the porous geometry can, however, not be derived. Studying how ${\textit{Bo}}_{{crit}}$ depends on, say, the minimum required distance between two obstacles would require completely new initial conditions for the trapped cluster geometries.

At this point, a comment on computational efficiency is appropriate, as we have repeatedly asserted that the quasistatic model is significantly faster than direct numerical simulations for predicting the critical Bond number. This claim may appear counterintuitive, given that the quasistatic model currently requires an initial condition that can only be obtained through a direct simulation. The key advantage, however, lies in the fact that only a single expensive simulation per cluster is needed to compute the critical Bond number as a function of the density ratio and, potentially, the contact angle. In contrast, without the quasistatic model, 5–10 direct simulations would be required for each combination of density ratio and contact angle. For a single pair of the density ratio and the contact angle, the speed-up is therefore roughly 5–10 times, while, if a detailed study across $(\rho _{{w}}/\rho _{nw}, \theta _0)$ -space was to be carried out, it is fair to say that employing the quasistatic model would lead to computational savings of several orders of magnitude.

In this connection, we note that the full numerical simulations based on phase-field modelling presented in §§ 4 and 5 are most likely not the most efficient way of initialising the quasistatic model. In the limit of small Bond number (that is, when surface tension is much more decisive for the shape of the cluster than the pressure gradient of the surrounding flow), the equilibrium configuration of the cluster surfaces could possibly be obtained directly, for example, by the Surface Evolver software of Brakke (Reference Brakke1992). The surrounding pressure field could subsequently be computed by solving the Stokes equation once for the entire domain employing the cluster boundaries found by the surface evolver. In the present work, we have, however, not tested this strategy and more research is needed to clarify how the initialisation of the quasistatic model can be accelerated.

4.1. DNS model

The phase-field model for two-phase flow, as formulated by Ding et al. (Reference Ding, Spelt and Shu2007), consists of a Navier–Stokes–Cahn–Hilliard system of differential equations supplemented with appropriate boundary conditions at solid interfaces. Letting $\boldsymbol{u}$ be the velocity of the fluids and $p$ be the pressure, the governing equations read

(4.1a) \begin{align} \rho (\phi ) \left [ \frac {\partial \boldsymbol{u}}{\partial t} + \big ( \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\big ) \boldsymbol{u} \right ] - \boldsymbol{\nabla }\boldsymbol{\cdot }\big ( 2 \mu (\phi ) \boldsymbol{D u} \big ) + \boldsymbol{\nabla }p & = g_\phi \boldsymbol{\nabla }\phi + \rho (\phi ) \boldsymbol{g}, \end{align}

(4.1b) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} & = 0, \end{align}

(4.1c) \begin{align} \frac {\partial \phi }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi & = \epsilon M_0 {\nabla} ^2 g_\phi , \end{align}

(4.1d) \begin{align} g_\phi & = \frac {3 \gamma }{2 \sqrt {2}} \Big ( \epsilon ^{-1} W'(\phi ) - \epsilon {\nabla} ^2 \phi \Big ). \end{align}

Here, $\rho (\phi )$ and $\mu (\phi )$ are linear interpolations of the mass densities and dynamic viscosities of the two fluids, respectively. For example,

(4.2) \begin{align} \rho (\phi ) \equiv \frac {1 + \phi }{2} \rho _{{w}} + \frac {1-\phi }{2}\rho _{{nw}}, \end{align}

with $\mu (\phi )$ defined completely analogously. Moreover, $\boldsymbol{Du} \equiv (\boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^\top )/2$ is the rate-of-strain tensor, $M_0$ is the phase field mobility and $W(\phi ) \equiv ( 1 - \phi ^2)^2 /4$ is the Ginzburg–Landau double-well potential. It should be noted that (4.1) is not strictly thermodynamically consistent; in contrast, the model by Abels et al. (Reference Abels, Garcke and Grün2012) obtains thermodynamic consistency by including a small term in (4.1a ) that makes the model respect a version of the second law of thermodynamics for finite $\epsilon \gt 0$ . However, for small $\epsilon$ , we expect little difference between the two formulations, as shown numerically by Aland & Voigt (Reference Aland and Voigt2012).

At solid boundaries, the following conditions apply (see e.g. Carlson, Bellani & Amberg Reference Carlson, Bellani and Amberg2012):

(4.3a) \begin{align} \boldsymbol{u} &= \boldsymbol{0}, \end{align}

(4.3b) \begin{align} \hat {\boldsymbol{n}} \boldsymbol{\cdot }\boldsymbol{\nabla }g_\phi &= {0}, \end{align}

(4.3c) \begin{align} \hat {\boldsymbol{n}} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi &= \frac {2\sqrt {2}}{3 \epsilon } \cos (\theta _0) f'_w (\phi ) , \end{align}

where $f_w(\phi ) \equiv (2 + 3\phi - \phi ^3)/4$ smoothly interpolates between 0 and 1 to enforce the difference in the solid–liquid interface energy between the two liquid phases, and $\hat {\boldsymbol{n}}$ is the normal vector of the solid boundary pointing into the solid. Periodic boundary conditions are applied to all fields on the horizontal and vertical sides of the domain.

We employ the open source finite element code Twoasis (available from the depository of Linga at https://github.com/gautelinga/twoasis) to solve the governing equations of the DNS model numerically. This code is based on the Oasis solver (Mortensen & Valen-Sendstad Reference Mortensen and Valen-Sendstad2015) built on the open source framework Fenics (see Logg, Mardal & Wells Reference Logg, Mardal and Wells2012). As such, it is capable of discretising the spatial dimensions of the governing equations using, e.g. Lagrange elements of any desired order. However, since Twoasis carries out the time-integration using a semi-implicit splitting scheme which is first-order accurate in the time step, we use first-order spatial elements to treat space and time equally. We generate doubly periodic triangular meshes for the simulations using the built-in meshing tool in the Bernaise solver (Linga, Bolet & Mathiesen Reference Linga, Bolet and Mathiesen2019). This tool uses the Python library MeshPy, which is based on the mesh generator Triangle written by Shewchuk (Reference Shewchuk1996). For further details about the code and its implementation, we refer the reader to the source code, which is freely available from Github at the URLs given earlier. The parameter values used in the DNS computations are specified in the next section.

4.2. Simulation protocol and computational parameters

We solve the governing equations (4.1) in a square domain of size $L \times L$ and start all our simulations with the velocity field, $\boldsymbol{u}$ , the pressure, $p$ , and the phase field chemical potential, $g_\phi$ , equal to zero everywhere in the domain. The values of $\phi$ are chosen such that the domain is initially filled with non-wetting fluid except for a circle of radius $r_{\textit{circ}}$ which contains wetting fluid. Naturally, then, $r_{\textit{circ}}$ becomes the length scale for the size of the trapped cluster. We emphasise that using $r_{\textit{circ}}$ as a measure for the cluster size disregards that a cluster, for example, may be longer than it is wide, as has been shown by Tallakstad et al. (Reference Tallakstad, Knudsen, Ramstad, Løvoll, Måløy, Toussaint and Flekkøy2009b ) for very large clusters (see the inset of their figure 3). For simplicity, we neglect this effect throughout this paper and note that a more detailed description of the stable cluster size would consider the cluster’s extension in both the $x$ - and $y$ -directions. The location of the circle’s centre is chosen randomly from a uniform distribution over the domain and, at later stages, we obtain statistical properties of the cluster stability for a given value of $r_{{circle}}/r$ by carrying out the simulation for an ensemble of random locations of the circle centre drawn from a uniform distribution over the domain. To treat clusters of different sizes similarly, we choose $L = 8 r_{\textit{circ}}$ ; in this way, the fraction of the domain which is blocked by the trapped cluster becomes independent of the value of the circle radius.

To enable comparisons, we use the same fluid properties as those used by Mathiesen et al. (Reference Mathiesen, Linga, Misztal, Renard and Le Borgne2023) in our simulations, that is, we let $\rho _{{w}}/\rho _{{nw}} = 1.25$ , $\mu _{{w}}/\mu _{{nw}} = 0.7125$ and $\theta _0 = \pi /3$ . These values correspond to a system filled with water and silicone oil. Moreover, we fix the absolute values of $r$ , $\rho _{{w}}$ , $\mu _{{w}}$ and $g$ , and let the width of the phase transition be $\epsilon /r = 0.05$ , the mesh size be $h/r = 0.05$ and the phase field mobility be $\epsilon M_0/(\rho _{{w}} (r/g)^{3/2}) = 10^{-11/2}$ . As is illustrated in figure 6(b), these choices ensure that the interface is resolved by 3–4 elements. This is enough to retain high accuracy in the simulations.

Figure 6.

(a) Initial and final states of a cluster for $r_{\textit{circ}}/r = 4$ and ${Bo = 0.025}$ . Over time, the wetting fluid is redistributed from the circular configuration to one where the contact angle is everywhere equal to $\theta _0$ . (b) Close-up on the diffuse interface, whose width is proportional to $\epsilon$ , in one of the cluster openings in the final state together with the employed triangular mesh. As can be seen, the interface is resolved by 3–4 elements.

Now, as was shown by Mathiesen et al. (see their figure 2a), when using the previously mentioned physical parameters, the average flow velocity, $\langle v \rangle$ , satisfies

(4.4) \begin{align} \langle v \rangle \lesssim \frac {\gamma }{\mu _{{w}}} \frac {{\textit{Bo}}}{100}. \end{align}

We therefore fix the velocity scale $v_0 \equiv \mu _{{w}}/(1000 \rho _{{w}} r)$ and choose the surface tension as $\gamma = 100 v_0 \mu _{{w}}/{\textit{Bo}}$ , for in that case, the average flow velocity will be close to $v_0$ , while the characteristic Reynolds number will be

(4.5) \begin{align} {\textit{Re}} \equiv \frac {\rho _{{w}} r v_0}{\mu _{{w}}} \lesssim 10^{-3}, \end{align}

ensuring that the flow is indeed creeping. We note that this choice implies that the Capillary number is ${\textit{Ca}} = O({\textit{Bo}}/100)$ , justifying the simplification of (3.14) to (3.15). In addition, the velocity scale $v_0$ induces the advection time scale $t_0 \equiv 2 r / v_0$ and this allows us to define a stopping criterion for the simulations. When the simulation is started and the cluster is allowed to evolve from the circular configuration, its geometry initially undergoes rapid changes before approaching a fixed shape, in which the only changes in the phase field are due to diffusion at the cluster interface. As such, it is reasonable to stop the simulations when the phase field has not changed substantially (defined through some threshold value) anywhere in the domain during the last time interval of length $t_0$ . In practice, before the start of a simulation, we draw a set $\{ \boldsymbol{r}_n \}$ of $1000 (L/r)^2$ uniformly distributed random points inside the domain and carry out the time integration up to $t = t_{\!f}$ , where $t_{\!f}$ is the smallest time at which the phase field satisfies the criterion

(4.6) \begin{align} \frac {1}{2} \max _{n} \Big ( \big | \phi (\boldsymbol{r}_n, t_{\!f}) - \phi (\boldsymbol{r}_n, t_{\!f} - t_0) \big | \Big ) \lt \frac {1}{100}. \end{align}

That is, we carry out the time integration until the absolute change in the phase field at any of the points in the set $\{ \boldsymbol{r}_n \}$ is less than $1/50$ during the latest time interval of length $t_0$ . Finally, we choose the time step based on a CFL-type condition. To be precise, we take the time step to be $\Delta t = C_r h/v_0$ , where $C_r \equiv 0.004/(r_{\textit{circ}}/r)$ is the Courant number. Based on experience, this is approximately the largest possible time step allowed by the numerical time integration that ensures stability. We have chosen the Courant number to be inversely proportional to the cluster size to be able to simulate lower Bond numbers for larger clusters. An example of a cluster’s initial and final shapes are shown in figure 6(a) for $r_{\textit{circ}}/r = 4$ , ${\textit{Bo}} = 0.025$ . The reader will notice that there appears to be slightly more wetting fluid in the final state than in the initial state. This is because the two fluids are not strictly incompressible in the phase field model. Instead, as a consequence of (4.1b ) and (4.1c ), the phase field model conserves the total phase field in the sense that

(4.7) \begin{align} \frac {\mathrm{d} }{\mathrm{d} t} \int _{\mathcal{D}} \phi \, \mathrm{d}^2 \boldsymbol{x} = 0, \end{align}

in which $\mathcal{D}$ is the fluid domain. Therefore, the phase field model may lead to changes in cluster size. As discussed by Yue, Zhou & Feng (Reference Yue, Zhou and Feng2007), this effect is particularly pronounced for large values of $\epsilon /r_{\textit{circ}}$ and $L / r_{\textit{circ}}$ . In addition, we have found that the surface–liquid interaction, enforced through (4.3c ), also influences the final size of the cluster. With the presently used parameters, the relative change in area is, however, small enough that we describe all cluster sizes in terms of their initial radii.

4.3. Comparing the quasistatic model to the DNS model

The purpose of the quasistatic model is to predict the critical Bond number at which a given cluster becomes unstable. The model will do so if it accurately predicts the evolution of the opening angles as a function of ${\textit{Bo}}$ , and we may therefore validate it by showing that its results for $\boldsymbol{\omega }({\textit{Bo}})$ are close to those of the DNS model. To do so, we first use the DNS model to compute the opening angles of two clusters for a discrete set of Bond numbers. One cluster starts out as a circle with $r_{\textit{circ}}/r=2$ , while the other is twice as large and has $r_{\textit{circ}}/r = 4$ . Since the interfaces in the DNS model are diffuse, we compute the opening angles from the intersection points between the obstacles and the curves along which $\phi = 0$ . When measured in radians, the uncertainty of the opening angles is therefore approximately given by the width of the phase transition, $\epsilon$ . Having solved the DNS model for a range of Bond numbers, we then initialise the quasistatic model from the opening angles produced by the DNS model for the smallest value of the Bond number and integrate (3.1) with respect to ${\textit{Bo}}$ using the classical fourth-order Runge–Kutta method with a step size of $\Delta {\textit{Bo}} = 0.0001$ . This value of $\Delta {\textit{Bo}}$ is small enough to guarantee that the integration errors are negligible relative to the model errors.

Figure 7.

(a,b) Phase field of a trapped cluster with $r_{\textit{circ}}/r = 2$ for (a) ${\textit{Bo}} = 0.05$ and (b) ${\textit{Bo}} = 0.2125$ . The dashed, black lines are the interfaces predicted by the quasistatic model. Since this model is initialised from the results of the DNS model for ${\textit{Bo}} = 0.05$ , the deviation between the models is zero at this Bond number. (c) Evolution of the opening angles, $\omega _n$ , normalised by their critical values given by (3.24) as a function of the Bond number. The full lines show the result of the quasistatic model, while the circles show the results of the DNS model. In the DNS computations, the uncertainty of $\omega _n$ is assumed to be $\pm \epsilon$ due to the diffusive interface, and this is indicated by the error bars. The numbering $n = 1, 2, \ldots , 5$ corresponds to the numbering of the cluster openings in panel (a).

The results for the cluster with $r_{\textit{circ}}/r = 2$ are shown in figure 7. The agreement between the two models can be seen to be very good, with essentially only a few data points of the DNS model for interfaces 2 and 3 close to the critical Bond number failing to agree with the quasistatic model within the uncertainty. In particular, the quasistatic model’s prediction of the evolution of $\omega _1$ , which is the opening angle causing the cluster to become unstable, is very close to the prediction of the DNS model. While the quasistatic model predicts the critical Bond number of this cluster to be ${\textit{Bo}}_{{crit}} = 0.222$ , we have checked that the DNS model predicts instability for any Bond number larger than $0.2375$ . As such, the quasistatic model’s prediction of the critical Bond number deviates by maximally 6.5 % from that of the DNS model, and we note that this deviation could be made smaller by carrying out additional simulations in the range $0.225 \leqslant {\textit{Bo}} \leqslant 0.2375$ .

For the cluster with $r_{\textit{circ}}/r = 4$ , the story is very much the same, as can be seen in figure 8. In general, the quasistatic model and the DNS model agree quite closely, and the few existing deviations only appear as the Bond number approaches its critical value, this time mainly for interfaces 6 and 8. For this particular cluster, the quasistatic model predicts that the critical Bond number is $0.1356$ and since the DNS model yields instability for any Bond number larger than $0.143$ , the deviation is maximally 4.5 %. Again, this deviation could be reduced by executing the DNS model for more values of the Bond number.

This concludes our validation of the quasistatic model. Although more cases can always be checked, we consider the two presented here to provide sufficient evidence that the quasistatic model can in fact accurately predict the critical Bond number of a given cluster.

Figure 8.

(a,b) Phase field of a trapped cluster with $r_{\textit{circ}}/r = 4$ for (a) ${\textit{Bo}} = 0.025$ and (b) ${\textit{Bo}} = 0.131$ . The dashed, black lines are the interfaces predicted by the quasistatic model. Since this model is initialised from the results of the DNS model for ${\textit{Bo}} = 0.025$ , the deviation between the models is zero at this Bond number. (c) Evolution of the opening angles, $\omega _n$ , normalised by their critical values given by (3.24) as a function of the Bond number. The full lines show the result of the quasistatic model, while the circles show the results of the DNS model. In the DNS computations, the uncertainty of $\omega _n$ is assumed to be $\pm \epsilon$ due to the diffusive interface, and this is indicated by the error bars. The numbering $n = 1, 2, \ldots , 11$ corresponds to the numbering of the cluster openings in panel (a).

Figure 8. Long description

Panel A: A phase field diagram with a trapped cluster for a Bond number of 0.025. The diagram features gray circles within a red and blue background, with dashed black lines indicating interfaces predicted by the quasistatic model. Panel B: Another phase field diagram with a trapped cluster for a Bond number of 0.131, similar in structure to Panel A. Panel C: A line graph showing the evolution of normalized opening angles as a function of the Bond number. The x-axis represents the Bond number, and the y-axis represents the normalized opening angles. The graph includes multiple colored lines and circles with error bars, indicating the results of the quasistatic model and the DNS model, respectively. The numbering corresponds to the cluster openings in Panel A.

5.1. Probability that a cluster is stable as a function of its size and the Bond number

To compute the probability that a cluster of a given size is stable for a given Bond number, we first construct initial conditions for the quasistatic model using the DNS model. For each of the cluster sizes $r_{\textit{circ}}/r = 2$ , 4, 6, 8 and 10, we employ the DNS model with ${\textit{Bo}} = (10 r_{\textit{circ}}/r)^{-1}$ to compute the stable configuration of 10 clusters starting at random positions in the domain. Each stable state is then used to initialise the quasistatic model, which is subsequently integrated up to its critical Bond number, ${\textit{Bo}}_{{crit}}$ . If we let $N_{\textit{stable}}(r_{\textit{circ}}/r, {\textit{Bo}})$ denote the number of clusters of size $r_{\textit{circ}}$ which have ${\textit{Bo}}_{{crit}} \gt {\textit{Bo}}$ , this enables us to compute the probability, $P$ , that a cluster of size $r_{\textit{circ}}/r$ is stable at Bond number ${\textit{Bo}}$ as

(5.1) \begin{align} P(r_{\textit{circ}}/r, {\textit{Bo}}) = \frac {N_{\textit{stable}}(r_{\textit{circ/r}}, {\textit{Bo}})}{10}. \end{align}

The choice of using only 10 trapped clusters for each size is a consequence of our available computing resources; using four Intel Xeon Gold 6148 CPUs having 20 cores each, it took us approximately four weeks to generate the initial conditions for the case $r_{\textit{circ}}/r = 10$ alone. Generating, for example, 100 initial conditions would hence be unfeasible without significantly larger computational resources. In contrast, the part of the computation executed using the quasistatic model takes less than a minute on a single CPU core even for the largest cluster size, the ‘bottleneck’ being the extraction of the opening angles from the DNS model. Hence, the computation of (5.1) would have been unfeasible even for a modest resolution in ${\textit{Bo}}$ without the quasistatic model.

Figure 9.

Probability that a cluster of size $r_{\textit{circ}}/r$ is stable at a given Bond number as a function of (a) ${\textit{Bo}}$ , (b) ${\textit{Bo}}^{1/2} \times (r_{\textit{circ}}/r)$ , (c) ${\textit{Bo}}^{4/5} \times (r_{\textit{circ}}/r)$ and (d) ${\textit{Bo}} \times (r_{\textit{circ}}/r)$ . The legend applies to all figures. According to MFA1 and MFA2, the results for $P(r_{\textit{circ}}/r, {Bo)}$ should collapse to a single curve when plotted as a function of ${\textit{Bo}} \times (r_{\textit{circ}}/r)$ and ${\textit{Bo}}^{1/2} \times (r_{\textit{circ}}/r)$ , respectively.

Now, our results for the probability of stability are shown in figure 9. When the Bond number is fixed, figure 9(a) illustrates the fact that large clusters are more prone to be unstable than small clusters, as one would expect. However, our results do not agree with the predictions of MFA1 and MFA2 described in § 1. According to MFA2, the probability (5.1) should be a function of ${\textit{Bo}}^{1/2} \times (r_{\textit{circ}}/r)$ only and the results for $P(r_{\textit{circ}}/r, {Bo)}$ should therefore collapse to a single curve when plotted against this parameter. As is clearly seen in figure 9(b), our results do not have this property. Likewise, as shown in figure 9(d), our results also do not follow the scaling predicted by MFA1, i.e. that the probability should be a function of the parameter ${\textit{Bo}} \times (r_{\textit{circ}}/r)$ only. The disagreement with MFA2 is, of course, particularly noteworthy, since we have used the same fluid properties as Mathiesen et al. (Reference Mathiesen, Linga, Misztal, Renard and Le Borgne2023), and because they found a very good agreement between MFA2 and their numerical simulations. It is, however, important to note that Mathiesen et al., simulated a dynamic flow in a very large domain where clusters could both split and merge continuously. Moreover, their domains contained a large number of both trapped and mobilised clusters, which they did not distinguish between when comparing to the analytical predictions. As such, the agreement found by Mathiesen et al. (Reference Mathiesen, Linga, Misztal, Renard and Le Borgne2023) may be induced by their dynamic setting.

By varying the exponent of ${\textit{Bo}}$ , we find that the best collapse of our data to a single curve is observed when it is plotted as a function of ${\textit{Bo}}^{4/5} \times (r_{\textit{circ}}/r)$ . As shown in the figure 9(c), the collapse is reasonably good for this exponent of ${\textit{Bo}}$ , which means that the size of the largest stable cluster should scale in proportion to $r/{\textit{Bo}}^{4/5}$ . However, this result only holds for the set of fluid properties that are presently used.

5.2. Shortcomings of the mean-field theories

Having shown that both MFA1 and MFA2 do not agree with the simulated size of the largest stable cluster, we now discuss some reasons for why this is the case.

5.2.1. Pressure field around a trapped cluster

In any mean-field theory, the pressure field surrounding the cluster is approximated through the average pressure gradient of the flow. For the present gravity-driven setting, this corresponds to the assumption that the difference between the pressure just outside of the $n$ th and $i$ th openings of the cluster is given by

(5.2) \begin{align} p_{\textit{out},n} - p_{\textit{out},i} = \rho _{{nw}} g (y_{\textit{out},n} - y_{\textit{out},i}). \end{align}

Unless the cluster is large, this equation will only provide a poor approximation for the pressure difference between its top and bottom. In essence, the explanation for this is that the pressure difference between two different openings is a random variable, whose standard deviation is independent of the difference between the two openings’ $y$ -coordinates. While (5.2) can in fact provide a reasonably good estimate of the pressure difference’s mean value, the mean-field approximation can only be considered accurate if the standard deviation is small relative to the mean value. Since this requires a large difference in $y$ -coordinates, the mean-field approximation can therefore only be accurate for large clusters.

To illustrate this point for the presently used porous geometry, figure 10 shows the pressure differences between the openings obtained from the initial conditions produced by the DNS model used for figure 9 as a function of their difference in $y$ -coordinate. From this, it is clear that the pressure difference is a random variable whose mean value is reasonably well described by (5.2), while its standard deviation is roughly $3 \rho _{{w}}gr$ independently of the difference in $y$ -coordinates. This means that if we want the standard deviation to be smaller than, say, $5\,\%$ of the mean-field prediction, we must have $(y_{\textit{out},n} - y_{\textit{out},1})/r \simeq 75$ . Thus, the pressure difference between the top and bottom of a particular cluster can only be expected to be within $5\,\%$ of the mean-field prediction if the cluster’s size is $r_{\textit{circ}}/r \gtrsim 32$ , which is indeed quite large. The inset in figure 10 shows the ratio between the standard deviation and the mean value of the pressure difference as a function of the difference in $y$ -coordinates. Clearly, this ratio diverges for small differences in $y$ -coordinates, and for none of the clusters considered in this study is the ratio smaller than $10\,\%$ .

Figure 10.

Pressure difference $\Delta p_{n1} \equiv p_{\textit{out},n} - p_{\textit{out},1}$ as a function of the difference in $y$ -coordinates between the evaluation points, $\Delta y_{n1} \equiv y_{\textit{out},n} - y_{\textit{out},1}$ . The green circles are obtained from the initial conditions produced by the DNS model using the outside evaluation point $\boldsymbol{r}_{\textit{out},n} \equiv \boldsymbol{r}_{\textit{int},n} - \delta \boldsymbol{\hat {n}}_n$ with $\delta = 3 \epsilon$ . The red, dashed line shows the mean-field approximation given by (5.2). The inset shows the standard deviation, $\sigma _{\Delta p}$ , of the green circles divided by the mean value, $\mu _{\Delta p}$ , of the pressure difference as a function of the difference in $y$ -coordinates. These quantities are computed from the green circles for each of the segments $m/2 \leqslant \Delta y/r \lt (m+1)/2$ for $m = 0, 1, 2, \ldots$ . Roughly, we have $\sigma _{\Delta p} \simeq 3 \rho _{{w}} g r$ independently of $\Delta y$ .

Figure 11.

The pressure field around a trapped cluster of size $r_{\textit{circ}}/r = 8$ .

Figure 11. Long description

A heat map showing the pressure field around a trapped cluster. The heat map uses a color gradient to represent the pressure values, with the color scale ranging from -35 to 10. The x-axis is labeled x/r and ranges from -15 to 15. The y-axis is labeled y/r and ranges from 0 to 20. The color gradient transitions from purple at the lowest values to blue, green, yellow, orange, and red at the highest values. The central region of the heat map shows higher pressure values, indicated by the red and orange colors, while the outer regions show lower pressure values, indicated by the blue and green colors. The heat map also includes circular markers distributed throughout the grid, with black dots at their centers.

The physical reason for the fixed standard deviation of the pressure difference is that, locally, the pressure field is much more affected by the geometry of the medium than by the average pressure gradient of the flow. This is illustrated by the trapped cluster depicted in figure 11, for which the mean-field assumption is clearly inaccurate. Instead, the pressure field remains approximately constant over, say, 5–10 cluster openings before changing abruptly because the surrounding fluid must go through a narrow passage. In porous geometries that either contain less narrow passages or exhibit a higher degree of connectivity, such as a close packing of spheres in three dimensions, there will be fewer of these abrupt changes and the standard deviation of the pressure differences between the cluster openings is therefore expected to be smaller. As a consequence, the mean-field approximation (5.2) is expected to be more accurate for these geometries than that used in the present work. For example, it may work well for all but the smallest clusters in the highly permeable media used by Tallakstad et al. (Reference Tallakstad, Løvoll, Knudsen, Ramstad, Flekkøy and Måløy2009a ).

5.2.2. Width of the breaking opening

In MFA1, it is assumed that a cluster breaks when the pressure difference (approximated through the average pressure gradient) across the cluster is larger than the maximum sustainable pressure difference of a typical opening. If this was true, one would expect that the width of the opening that actually breaks on average is the same as the average width of a random cluster opening. As illustrated in figure 12, the average width of the breaking opening increases from $\langle d_{\textit{break}}/r \rangle \simeq 1.0$ for small clusters to $\langle d_{\textit{break}}/r \rangle \simeq 1.7$ for large clusters. In contrast, the average width of a random opening is nearly independent of the cluster size and has $\langle d_{\textit{random}}/r \rangle \simeq 0.8$ . As such, this assumption made in MFA1 is not generally accurate.

Interestingly, it is assumed in MFA2 that a cluster breaks if its widest opening cannot sustain the pressure difference across the cluster. Considering the results shown in figure 12, the width of the breaking opening is on average clearly smaller than the average width of the widest opening. However, the breaking opening is most often located on the lower half of the cluster and the average width of the widest opening on the lower half of the cluster is typically quite close to the average width of the breaking opening. In that sense, the assumption made in MFA2 can be said to be satisfied.

Figure 12.

Average width of the breaking opening, the widest opening, the widest opening on the lower half of the cluster as well as a random opening as functions of the cluster size $r_{\textit{circ}}/r$ . The uncertainty indicated by the error bars is the estimated standard deviation divided by the square root of the number of samples.

Figure 13.

Probability density function of the non-dimensional maximum sustainable pressure, $\eta$ , for a random cluster opening given by (5.4) for the contact angle $\theta _0 = \pi /3$ . The PDF of $d/r$ in this expression is sampled from the edge lengths of a Delaunay triangulation of the porous geometry, and the uncertainty induced by the finite number of edges is what causes the small, local oscillations of the curve. The cutoff indicated by the dashed, black line is a consequence of the requirement that the distance between any two obstacles, $d$ , must be greater than the minimum distance $d_{\textit{min}} = r/10$ . The dashed, red line is the function $(\eta \times {\textit{Bo}})^{-7/4}$ and shows that $p_\eta$ asymptotically (but prior to the cutoff) follows a power law.

5.2.3. Probability density function of the maximum sustainable pressure difference of a random opening

In § 1, we remarked that the exponential distribution for the maximum sustainable pressure of a random opening assumed in MFA2 is not supported by theory. Here, we go one step further and argue that the assumption is in fact incorrect. To that end, we start by defining the maximum non-dimensional sustainable pressure difference $\eta \equiv \text{max} (\Delta p /(\rho _{{w}} g r))$ . Inserting (3.24) into (3.23) leads to

(5.3) \begin{align} \eta = \frac {2}{{\textit{Bo}}} \frac {1}{\big ( (d/r + 2)^2 - 4 \sin ^2(\theta _0) \big )^{1/2} - 2\cos (\theta _0)} \end{align}

for an opening of width $d$ . Since the opening width is a random variable, $\eta$ is likewise a random variable, and one can show that its probability density function (PDF) is given by

(5.4) \begin{align} p_\eta (\eta ) = \frac {{\textit{Bo}} \big ( 2 + 2\cos (\theta _0) \eta {\textit{Bo}} \big )}{(\eta {\textit{Bo}})^2 \big ( 1 + 2\cos (\theta _0) \eta {\textit{Bo}} + (\eta {\textit{Bo}})^2\big )^{1/2}} \times p_{d/r} \Bigg ( \frac {2}{\eta {\textit{Bo}}} \big ( 1 + 2\cos (\theta _0) \eta {\textit{Bo}} + (\eta {\textit{Bo}})^2 \big )^{1/2} - 2 \Bigg ), \end{align}

in which $p_{d/r}$ is the PDF of $d/r$ . To the best of our knowledge, $p_{d/r}$ is not known analytically and the amount of data that we have available for each value of $r_{\textit{circ}}/r$ does not allow us to sample it accurately. We therefore resort to an approximation and use the PDF of the edge length (minus $2r$ ) from Delaunay triangulations of the obstacle centres of a large number of porous geometries. Therefore, one should note that $p_{d/r}$ is expected to be smaller for small values of $d/r$ when based on the Delaunay triangulation than on actual trapped clusters. This is so because a trapped cluster tries to minimise its surface length to reach a local energy minimum and because the Delaunay algorithm tries to avoid short edges to maximise the minimum angle of the triangulation. As a consequence, $p_\eta$ is expected to have a heavier tail (that is, to decay less rapidly) when based on actual trapped cluster openings than on Delaunay edges. Our result for the PDF of $\eta$ is shown in figure 13 for the contact angle $\theta _0 = \pi /3$ . From the figure, it is clear that $p_\eta$ falls off as ${\textit{Bo}}^{-3/4} \times \eta ^{-7/4}$ asymptotically when based on Delaunay edges, and when combined with the previous remark, it is safe to conclude that $p_\eta$ does not fall off exponentially as a function of $\eta$ if it (more accurately) had been computed from actual trapped cluster openings. As such, the assumption in MFA2 made about the form of $p_\eta$ cannot be correct. It is interesting to note that if MFA2 is carried out with a distribution of the form $p_\eta (\eta ) \propto {\textit{Bo}}^{\beta +1} \times \eta ^{\beta }$ (which is implied by (5.4) for large values of ${\textit{Bo}} \times \eta$ if $p_{d/r}$ follows a power law distribution), one finds that the largest stable cluster size is very well approximated by $r/{\textit{Bo}}$ for the relevant range of Bond numbers regardless of the value of $\beta$ , as long as $\beta \lt -1$ . Hence, if the maximum sustainable pressure difference of a random opening follows a power law distribution, MFA2 in fact gives the same result as MFA1.

Figure 14.

Value of $\alpha$ as a function of the density ratio, $\rho _{{w}}/\rho _{{nw}}$ that leads to the best possible data collapse of the probability that a cluster of size $r_{\textit{circ}}/r$ is stable when plotted as a function ${\textit{Bo}}^\alpha \times (r_{\textit{circ}}/r)$ .

5.3. The $r/{\textit{Bo}}^{4/5}$ - scaling of the largest stable cluster size is not a general result

The value $4/5$ for the exponent of ${\textit{Bo}}$ found to give the best data collapse in figure 9 is implied by the particular values of the fluid properties. It is not a general result valid over a large part of the system’s parameter space.

To back up this claim, we consider the exponent’s dependency on the density ratio. To compute the critical Bond number for an arbitrary density ratio, we have first integrated the states computed by the DNS model for $\rho _{{w}}/\rho _{{nw}} = 1.25$ used to produce figure 9 up to different values of $\rho _{{w}}/\rho _{{nw}}$ using (3.25) before integrating the resulting states up to the critical Bond numbers using (3.1). Our results for the exponent, $\alpha$ , which gives the best collapse of the probability that a cluster of size $r_{\textit{circ}}/r$ is stable when plotted against the parameter ${\textit{Bo}}^{\alpha } \times (r_{\textit{circ}}/r)$ , are shown in figure 14 in the range $1.25 \leqslant \rho _{{w}}/\rho _{{nw}} \leqslant 10.0$ . From this, it is evident that the exponent is a non-trivial function of the density ratio. While we do not know the reason for the local minimum of $\alpha$ around $\rho _{{w}}/\rho _{{nw}} \simeq 5$ , it is relatively simple to understand why $\alpha$ approaches a value close to $1$ as the density ratio becomes large. In the limit of large density ratio, $p_{\textit{in},n} - p_{\textit{in},i}$ is much larger in magnitude than $p_{\textit{out},n} - p_{\textit{out},i}$ for any two openings $n$ and $i$ . As an approximation, one may therefore treat the pressure field outside the cluster as being constant in this limit. This means that if two openings have the same width, the one of the two located the farthest upstream will break first. For that reason, it is somewhat reasonable to assume that, among all openings, it is the one lying the most upstream that breaks first. Since the width of the most upstream opening is distributed just like any other opening, it also has the same average properties and the limit of large density ratio therefore approximately corresponds to the case analysed in MFA1. One therefore ends up with $\alpha$ being close to $1$ when the density ratio is large.

It is interesting to note that the prediction of MFA1 is expected to hold when $\rho _{{w}}/\rho _{{nw}}$ is very small and the medium is highly permeable. The explanation for this is as follows: when the trapped fluid is much lighter than the surrounding fluid, the internal pressure field becomes negligible. Moreover, if the medium is highly permeable, the surrounding pressure field is well approximated by the average pressure gradient of the flow. As a consequence, if two openings have the same width, the one of the two located the farthest upstream will break first. The argument used earlier then implies that $\alpha$ will be close to $1$ . It is tempting to explain the good match between MFA1 and the experiments performed by Tallakstad et al. (Reference Tallakstad, Løvoll, Knudsen, Ramstad, Flekkøy and Måløy2009a ) using this argument because the density ratio in the experiments is indeed very small and the medium is highly permeable. The explanation is, however, questionable since the experiments were performed in a horizontal cell, meaning that gravitational effects were absent.

Similar quantitative analyses as that mentioned could, of course, be carried out for $\alpha$ ’s dependency on other parameters such as the contact angle, $\theta _0$ , the minimum allowed distance between two obstacles, $d_{\textit{min}}$ , and other geometrical variables. While we have not done so, it seems unlikely that $\alpha$ should be independent of these variables. For example, as we have already argued, the minimum distance is a key parameter when it comes to the adequacy of pressure-field description in the mean-field arguments. When this distance is large, the optimal exponent will likely approach the prediction of the mean-field theory, while it remains unknown what happens in the limit where this distance is small. Classifying cluster stability in terms of a scaling relation of the form $r/{\textit{Bo}}^{\alpha }$ therefore requires that $\alpha$ is allowed to depend on several variables, leaving the efficiency of this representation questionable.