Interfacial dynamics and pinch-off singularities for axially symmetric Darcy flow

TL;DR

This paper investigates the development of pinch-off singularities in axially symmetric Darcy flow, revealing new behaviors not seen in 2D Hele-Shaw flow, supported by a novel level set numerical scheme and similarity analysis.

Contribution

It introduces a novel numerical method and demonstrates the occurrence of pinch-off singularities with a power-law behavior in axially symmetric porous media flow, extending understanding beyond 2D models.

Findings

Pinch-off singularities occur in 3D axially symmetric Darcy flow.

Minimum radius follows a power law with exponent 1/3 near pinch-off.

Differences between 3D and 2D flow behaviors are highlighted.

Abstract

We study a model for the evolution of an axially symmetric bubble of inviscid fluid in a homogeneous porous medium otherwise saturated with a viscous fluid. The model is a moving boundary problem that is a higher-dimensional analogue of Hele-Shaw flow. Here we are concerned with the development of pinch-off singularities characterised by a blow-up of the interface curvature and the bubble subsequently breaking up into two; these singularities do not occur in the corresponding two-dimensional Hele-Shaw problem. By applying a novel numerical scheme based on the level set method, we show that solutions to our problem can undergo pinch-off in various geometries. A similarity analysis suggests that the minimum radius behaves as a power law in time with exponent $\alpha = 1/3$ just before and after pinch-off has occurred, regardless of the initial conditions; our numerical results support this prediction. Further, we apply our numerical scheme to simulate the time-dependent development and translation of axially symmetric Saffman-Taylor fingers and Taylor-Saffman bubbles in a cylindrical tube, highlighting key similarities and differences with the well-studied two-dimensional cases.

Equations

$$\mbox{\boldmath$q$}=-\frac{K}{\mu}\mbox{\boldmath$\nabla$}p,$$

$$\displaystyle\mbox{\boldmath$v$}=\frac{\mbox{\boldmath$q$}}{\varphi},$$

$$p=p_B(t)-\gamma \kappa \text{mbox}on\partial \Omega ,$$

$$\begin{split}\mbox{\boldmath$x$}=L\hat{\mbox{\boldmath$x$}},\quad t=\frac{L^{3}\varphi}{Q}\hat{t},\quad\mbox{\boldmath$q$}=\frac{Q}{L^{2}}\hat{\mbox{\boldmath$q$}},\\ \mbox{\boldmath$v$}=\frac{Q}{\varphi L^{2}}\hat{\mbox{\boldmath$v$}},\quad p=p_{B}-\frac{\mu Q}{KL}\phi,\end{split}$$

$${\nabla }^2\varphi$$

$$v_n$$

$$\varphi$$

$$\frac{\partial \varphi }{\partial r}$$

$$\frac{1}{\beta }\frac{\partial \varphi }{\partial t}={\nabla }^2\varphi \text{in }{\mathbb{R}}^3\backslash \Omega (t),$$

$$\varphi \sim {\varphi }_{\infty }\text{as }r\to \infty$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial \varphi }{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial \varphi }{\partial \theta }\right)=0$$

$$s_t={\varphi }_r-\frac{1}{s^2}{\varphi }_{\theta }{s}_{\theta }$$

$$\varphi =\sigma \frac{3s{s}_{\theta }^2-\cot \theta s_{\theta }^3-s^2\left(s_{\theta \theta }+s_{\theta }\cot \theta \right)+2s^3}{s(s^2+s_{\theta }^2{)}^{3/2}}$$

$$\frac{\partial \varphi }{\partial r}\sim -\frac{1}{r^2}$$

$$\displaystyle\psi_{t}+F|\mbox{\boldmath$\nabla$}\psi|=0.$$

$$\displaystyle F=\frac{\mbox{\boldmath$\nabla$}\phi\cdot\mbox{\boldmath$\nabla$}\psi}{|\mbox{\boldmath$\nabla$}\psi|}\qquad\mbox{\boldmath$x$}\in\mathbb{R}^{3}\backslash\Omega.$$

$$\displaystyle\nabla^{4}F=0\qquad\mbox{\boldmath$x$}\in\Omega(t).$$

$$\frac{{\partial }^2\varphi }{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial \varphi }{\partial \rho }+\frac{{\partial }^2\varphi }{\partial z^2}=0,$$

$$\frac{\partial h}{\partial t}+\frac{\partial \varphi }{\partial z}\frac{\partial h}{\partial z}=\frac{\partial \varphi }{\partial \rho }$$

$$\varphi =\sigma \left(\frac{1}{h\sqrt{1+h_z^2}}-\frac{h_{zz}}{(1+h_z^2{)}^{3/2}}\right)$$

$$\varphi (\rho ,z,t)=(t_0-t{)}^{\gamma }\Phi (\xi ,\zeta ),h(z,t)=(t_0-t{)}^{\alpha }f(\zeta ),$$

$$\xi =\frac{\rho }{(t_0-t{)}^{\alpha }},\zeta =\frac{z-z_0}{(t_0-t{)}^{\beta }}.$$

$$\frac{{\partial }^2\Phi }{\partial {\xi }^2}+\frac{1}{\xi }\frac{\partial \Phi }{\partial \xi }+\frac{{\partial }^2\Phi }{\partial {\zeta }^2}=0,$$

$$-\frac{1}{3}f+\frac{1}{3}\zeta \frac{\mathrm{d}f}{\mathrm{d}\zeta }+\frac{\partial \Phi }{\partial \zeta }\frac{\mathrm{d}f}{\mathrm{d}\zeta }=\frac{\partial \Phi }{\partial \xi }$$

$$\Phi =\sigma \left(\frac{1}{f\sqrt{1+f^{\prime 2}}}-\frac{f^{\prime \prime }}{(1+f^{\prime 2}{)}^{3/2}}\right)$$

$$\alpha =1/3,\beta =1/3,\gamma =-1/3.$$

$$f^{\prime }\to f/\zeta \zeta \to \pm \infty .$$

$$\begin{array}{rl}p-p_B& =-{\left(\frac{\mu \phi {\gamma }^2}{K}\right)}^{1/3}(t_0-t{)}^{-1/3}\Phi ,\\ (\rho ,z-z_0,h)& ={\left(\frac{\gamma K}{\mu \phi }\right)}^{1/3}(\xi ,\zeta ,f).\end{array}$$

$$\varphi (\rho ,z,t)={\int }_{z_0-L}^{z_0+L}\frac{C(\hat{z},t)}{\sqrt{(\hat{z}-z{)}^2+{\rho }^2}}\text{d}\hat{z},$$

$$\begin{array}{r}\frac{\partial h}{\partial t}-\frac{\partial h}{\partial z}{\int }_{z_0-L}^{z_0+L}\frac{(z-\hat{z})C(z^{\prime },t)}{((\hat{z}-z{)}^2+h(z,t{)}^2{)}^{3/2}}\text{d}\hat{z}=\\ -h{\int }_{z_0-L}^{z_0+L}\frac{C(\hat{z},t)}{((\hat{z}-z{)}^2+h(z,t{)}^2{)}^{3/2}}\text{d}\hat{z},\end{array}$$

$$\begin{array}{r}{\int }_{z_0-L}^{z_0+L}\frac{C(\hat{z},t)}{\sqrt{(\hat{z}-z{)}^2+h(z,t{)}^2}}\text{d}\hat{z}=\\ \frac{1}{h\sqrt{1+h_z^2}}-\frac{h_{zz}}{(1+h_z^2{)}^{3/2}}.\end{array}$$

$$\frac{1}{3}(-f+\zeta f^{\prime })={\int }_{-\infty }^{\infty }\frac{(f^{\prime }(\zeta )(\zeta -\hat{\zeta })-f(\zeta ))D(\hat{\zeta })}{((\hat{\zeta }-\zeta {)}^2+f(\zeta {)}^2{)}^{3/2}}\text{d}\hat{\zeta },$$

$$\begin{array}{r}{\int }_{-\infty }^{\infty }\frac{D(\hat{\zeta })}{\sqrt{(\hat{\zeta }-\zeta {)}^2+f(\zeta {)}^2}}\text{d}\hat{\zeta }=\\ \frac{1}{f\sqrt{1+f^{\prime 2}}}-\frac{f^{\prime \prime }}{(1+f^{\prime 2}{)}^{3/2}}.\end{array}$$

$$\frac{1}{3}(-f+\zeta f^{\prime })=-\frac{D}{f}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}\hat{\eta }}{(1+{\hat{\eta }}^2{)}^{3/2}},$$

$$\begin{array}{r}\frac{1}{6}{\int }_{-\infty }^{\infty }\frac{f(\hat{\zeta })(f(\hat{\zeta })-\hat{\zeta }{f}^{\prime }(\hat{\zeta }))}{\sqrt{(\hat{\zeta }-\zeta {)}^2+f(\zeta {)}^2}}\text{d}\hat{\zeta }=\\ \frac{1}{f\sqrt{1+f^{\prime 2}}}-\frac{f^{\prime \prime }}{(1+f^{\prime 2}{)}^{3/2}}.\end{array}$$

Full text

PACS numbers

47.56.$+$r, 47.15.gp, 47.55.df, 47.55.N-, 47.55.Mh

††preprint: APS/123-QED

I Introduction

Studies of interfacial flows in porous media are motivated by various important applications such as oil recovery and salt-water intrusion [1]. There is particular interest in the physics community in studying immiscible flows that involve a fluid of lower viscosity displacing a fluid of higher viscosity since, in this scenario, the interface is unstable and prone to viscous fingering-type instabilities [2]. This paper provides a theoretical and numerical analysis of a one-phase model for flow in a homogeneous porous medium where motion of the more viscous fluid is governed by Darcy’s law

[TABLE]

while any pressure drop due to the less viscous (or inviscid) fluid’s motion is ignored. Here $q$ is the Darcy velocity, which can be thought of as a locally averaged velocity, while $p$ is the fluid pressure. The two constants $K$ and $\mu$ are the porous medium’s permeability and the fluid’s viscosity, respectively. The velocity of the fluid is related to the Darcy velocity via

[TABLE]

where $\varphi$ is the porosity. We shall be concerned with flows for which there is a region of inviscid fluid occupying the domain $\Omega(t)$ (which we refer to as a bubble), so that the boundary $\partial\Omega(t)$ is a moving interface to be determined as part of the problem.

As part of the problem specification, an interfacial condition on the pressure is required. We apply the macroscopic surface tension condition

[TABLE]

where $p_{B}$ is the pressure of the inviscid bubble, $\gamma$ is a dimensional surface tension and $\kappa$ is the mean curvature of the interface. While such a condition is used extensively in the context of two-dimensional Hele-Shaw models (which we shall discuss at length below), some justification of (3) is required. In real porous media, capillarity acts at the pore scale; for invasion by a non-wetting fluid, this leads to fractal-like fingering patterns on the scale of the individual pores, for which a macroscale, homogenised model is inadequate. On the other hand, there is significant experimental and theoretical evidence that for an advancing wetting fluid, the system exhibits an effective surface tension effect on a scale much larger that the pore size (see [3, 4, 5] and references therein). Macroscopic surface tension (3) has been used in studies of interfacial porous media flow since Chuoke et al. [6], and more recently by many others [7, 8, 9, 10, 11, 12, 13, 14]. We note that due to its origin, the boundary condition (3) may be more appropriate when an interface is advancing, but not receding, or vice versa. However, due to its simplicity and ubiquity in other studies, we will use it here in both cases. In the context of the fluid-fluid interface undergoing pinch-off, assuming the viscous fluid is the wetting fluid, the viscous fluid is advancing near the neck and so the inclusion of surface tension in (3) is reasonable.

We are particularly interested in scenarios for which the bubble $\Omega(t)$ develops a thin neck that subsequently pinches off and breaks into two parts. This change in topology involves a singularity in curvature of $\partial\Omega$, the type of which has been studied at length both mathematically and experimentally for a range of phenomena in fluid mechanics, including the break-up of inviscid bubbles [15, 16, 17]. Typically these studies are centred around seeking a self-similar form describing the collapse of the neck thickness by the power law $h_{\min}\propto(t_{0}-t)^{\alpha}$, where $t_{0}$ is the time at which pinch off occurs, and the similarity exponent $\alpha$ is a measure of the rate at which the velocity of the interface blows up. We refer the reader to [18, 19, 20] for a comprehensive overview of the analysis of these types of finite-time singularities in partial differential equations.

A further motivation for our study is to compare and contrast three-dimensional results (or at least axisymmetric results) for immiscible flow in a porous medium with the well-studied two-dimensional problem of flow through a Hele-Shaw cell. The analogy between these two problems arises because flow of a viscous fluid through a Hele-Shaw cell is also governed by Darcy’s law (1), except that the parameter $K$ in (1) is replaced by $b^{2}/12$, where $b$ is the small separation distance between the two parallel Hele-Shaw plates. To make direct comparisons, we consider (axisymmetric versions of) flow geometries that have received significant interest in the Hele-Shaw case. These include bubble contraction due to inviscid fluid being extracted from a point [21, 22, 23], displacement of viscous fluid in a channel geometry leading to Saffman-Taylor fingers [24, 25], propagation of finite (Taylor-Saffman) bubbles along a channel [26, 27], and injection of inviscid fluid at a point leading to tip-splitting and viscous fingering patterns in a radial geometry [28, 29, 30].

An additional observation to make about (3) is that our porous media flow model is also relevant for Stefan problems which describe certain melting or freezing phenomena; in that context, the boundary condition (3) is appropriate (here $p$ would need to be interpreted as temperature) as it describes the Gibbs-Thomson condition that relates melting/freezing temperature to the curvature of the solid-melt interface [31, 32, 33, 34, 35].

A further theoretical motivation for our work is to study the effect that the extra component of curvature has on the flow when we move from two dimensions to an axially symmetric geometry. To take an example, we note that while Hele-Shaw bubbles can undergo pinch-off [23, 22], this change in topology does not involve the same singularities in curvature as it does for three-dimensional flows (in two dimensions, the curvature of the interface does not blow up in the limit that pinch-off is approached); as such, a number of fundamental features of our examples will not be simple extensions of the Hele-Shaw cases. From a mathematical perspective, these issues are reminiscent of the differences between curve shortening flow and mean curvature flow in differential geometry [36, 37]. Given that Hele-Shaw flow may be interpreted as a nonlocal version of curve shortening flow in the plane [38, 21, 39], a final reason to use (3) is that our three-dimensional moving boundary problem can be considered a nonlocal version of the mean curvature flow. Again, in this context we are interested in studying the effect of surface tension in (3) on the development of axially symmetric curvature singularities.

We begin our article in Sec. II by summarising our model and explaining the connection to Stefan problems. The four flow geometries we consider are illustrated in Fig. 1. In Sec. III we treat the first of these geometries which involves a contracting bubble surrounded by an infinite body of fluid, as shown in Fig. 1$(a)$. This radial configuration is used extensively in the Hele-Shaw literature. Our numerical scheme, based on the level set method, is summarised for this radial geometry. The contracting bubble can eventually break up into two if the initial condition is sufficiently nonconvex. Following the methodology of Eggers et al. [15] normally associated with slender-body theory, we derive a similarity solution to describe the near pinch-off behaviour, predicting a similarity exponent $\alpha=1/3$. Our numerical solutions at times just before and after pinch-off follow the same similarity scaling and match well with the actual similarity solution in the former case. In Sec. IV we treat the second and third geometries in Fig. 1. These involve flow in a cylindrical channel, a set-up that is deliberately chosen as a higher-dimensional analogue of flow in a Hele-Shaw channel. We apply our (suitably adapted) numerical scheme and explore sufficiently nonconvex initial conditions that give rise to pinch-off singularities and convex initial conditions that evolve smoothly to travelling wave solutions. Such travelling wave solutions are axially symmetric analogues of the two-dimensional Saffman-Taylor fingers and Taylor-Saffman bubbles and thus our time-dependent study complements that of Levine & Tu [13] (see also Brener [10]), which was for steadily moving (axially symmetric) Saffman-Taylor fingers only. Finally, our results are summarised and discussed in Sec. V, where we briefly touch on the case of an expanding bubble [9, 12] (see Fig. 1$(d)$) which, provided a single finger is orientated along the axis of symmetry, can also lead to the same type of pinch-off as in the other geometries just mentioned.

II Model Formulation

II.1 Flow through a homogeneous porous medium

Our main physical motivation is to model the withdrawal of a bubble of inviscid fluid from a homogeneous porous medium that is otherwise saturated with viscous fluid. By assuming the viscous fluid is incompressible, \mbox{\boldmath\nabla}\cdot\mbox{\boldmathq}=0, then from (1) we have Laplace’s equation $\nabla^{2}p=0$ in the flow domain. We scale space, time, the Darcy velocity, the velocity of the viscous fluid, and pressure with

[TABLE]

where $\phi$ is the velocity potential satisfying \mbox{\boldmathq}=\mbox{\boldmath\nabla}\phi, $L$ is a representative length scale (related to the initial bubble shape or, for a more general configuration, the geometry of the flow) and $Q$ is the flow rate at which the inviscid fluid is injected/withdrawn. Retaining our original variable names, the resulting one-phase model for a bubble in a porous media is

[TABLE]

where $\sigma=\gamma K/(\mu Q)$ is the non-dimensional surface tension parameter and the positive and negative signs in (5d) correspond to the injection and withdrawal of the inviscid bubble, respectively. The boundary conditions on the bubble boundary $\partial\Omega(t)$ (5b) and (5c) are often referred to as kinematic and dynamic conditions, respectively. The kinematic condition (5b) simply states that the normal velocity of the boundary is the normal velocity of the fluid. The dynamic condition (5c) relates the pressure to the surface tension $\sigma$ via the mean curvature $\kappa$. This equation, the dimensionless version of (3), is discussed in the Introduction. From a geometric perspective, the surface tension acts to smooth out regions of high curvature on the interface. For the above mathematical model to make sense in the context of contracting bubbles, the inviscid fluid must be extracted from the point to which the interface ultimately evolves (or from multiple points in the case in which the initial bubble breaks up into two or more satellite bubbles).

For the special case of zero surface tension, $\sigma=0$, there has been a number of studies on contracting bubbles, focussing on the shape of the interface in the limit that it contracts to a point [40, 41, 42, 43], or whether the bubble breaks up into two [43] (in two dimensions, see [23, 44, 22, 45]). However, there has been relatively little theoretical study concerned with solutions to (5a)-(5d) that undergo changes in topology when the effects of surface tension are included. For expanding bubbles, with the positive sign in (5d), or for the problem involving displacement of viscous fluid in a cylindrical tube, details of linear, weakly nonlinear and travelling wave analysis has been performed [13, 12, 9]. However, for these geometries, no fully nonlinear time-dependent solutions to (5a)-(5d) (with appropriate changes due to relevant configuration) have been previously reported.

II.2 Application to Stefan problems

The other physical motivation for (5a)-(5d) is to model the melting/freezing of a crystal dendrite, where now $\partial\Omega$ represents a solid-melt interface. Typically, melting/freezing problems are modelled mathematically as a Stefan problem, where $\phi$ represents temperature, and the field equation is the linear heat equation

[TABLE]

where $\beta$ is the Stefan number, which is the ratio of latent heat to specific heat. In the regime for which the latent heat is significantly larger than the specific heat, it follows that $\beta\gg 1$ so that (6) reduces to (5a) [46]. In this context, the far-field condition (5d) represents the heat flux being prescribed at infinity. Another, perhaps more suitable, boundary condition in the context of Stefan problems that the temperature in the far field is prescribed leading to

[TABLE]

[43]. As we demonstrate a number of times throughout the paper, the dynamics of the interface near pinch-off are dependent only upon (5a)-(5c), and therefore independent of the choice of far-field boundary condition.

Pinch-off of crystals undergoing phase change is relevant, for example, for experiments conducted by Ishiguro et al. [47]. In this experiment, 3He crystals were grown in a 4He bath, and the authors argue the effects of gravity and surface tension cause the 3He crystal to pinch-off, as illustrated in Fig. 2. It was shown that the neck of the crystal both moves towards and recoils from pinch-off with a similarity exponent of $\alpha=1/3$. We return to this point in Sec. V.

III Contracting bubble geometry

III.1 Numerical scheme

In this section, we are concerned with the behaviour of a contracting bubble that undergoes pinch-off. Restricting ourselves to axially symmetric geometries in spherical co-ordinates such that $\phi=\phi(r,\theta,t)$ and with the bubble interface $\partial\Omega$ denoted by $r=s(\theta,t)$, our model (5a)-(5d) is

[TABLE]

The simulation of interfaces undergoing a topological change can be a computationally demanding task. For irrotational, inviscid flow, perhaps the most popular choice of numerical scheme is the boundary integral method (BIM) [48, 49, 50, 51]. While the BIM has shown to produce accurate simulations of interfaces approaching pinch-off, it can only be used for problems where the Green’s functions can be computed. Furthermore, on its own the BIM cannot easily handle changes in topology. Success has been achieved using ‘front capturing methods’, which are generally not as accurate but have greater flexibility [52, 53, 54]. The numerical scheme we use in this article implements the level set method [55], a numerical framework that describes the evolution of moving interfaces by representing them as the zero level set of a higher dimensional hyper-surface. Our scheme was presented in [43], and used to study solutions to (5a)-(5c) and (7) that undergo pinch-off for the special case of zero surface tension, $\sigma=0$. In that study it was shown that our numerical scheme is capable of accurately describing the behaviour of interfaces that contract to either one or multiple points. In this section, we give a brief description of the scheme, and refer the reader to [43] for further details.

To implement the level set method, we construct a signed distance function, \psi(\mbox{\boldmathx},t), where $\psi<0$ in the bubble region and $\psi>0$ in the viscous fluid region. If our bubble interface has a normal speed $V_{n}$, then we wish to construct a continuous function $F$ which satisfies $F=V_{n}$ on the interface. Thus $\psi$, and in turn the interface, evolves according to

[TABLE]

In the context of (5a)-(5d)

[TABLE]

The spatial derivatives in (9) are approximated using a second order essentially non-oscillatory scheme, and integrated in time using third order Runge-Kutta. To ensure the numerical solution accurately captures the dynamics of the interface both as it approaches and recoils from pinch-off, we choose a time-step size $\Delta t=0.05\times\Delta x/\max|F|$. Simulations are performed on the computational domain $0\leq r\leq 1.5$ and $0\leq\theta\leq\pi$ using $300\times 630$ equally spaced nodes.

We solve for $\phi$ using a similar procedure as described in [56]. For nodes away from the interface, the derivatives in (8a) are evaluated using a standard five-point finite difference stencil. For nodes adjacent to the interface, the stencil is adjusted by imposing a ghost node on the interface, whose location is determined by computing points where $\psi=0$. The value of this ghost node comes from the dynamic boundary condition (8c), where the curvature term is computed using \kappa=\mbox{\boldmath\nabla}\cdot\left(\mbox{\boldmath\nabla}\psi/|\mbox{\boldmath\nabla}\psi|\right). We incorporate the far-field boundary condition (8d) into the stencil by implementing a Dirichlet-to-Neumann mapping [57]. While solving for $p$ allows us to compute $F$ where \mbox{\boldmathx}\in\mathbb{R}^{3}\backslash\Omega, to solve (9) we require an expression for $F$ over the entire computational domain. It was proposed in [58] that $F$ can be extended into \mbox{\boldmathx}\in\Omega by solving the biharmonic equation

[TABLE]

Computing the solution to (11) does not require the location of the interface explicitly. Instead the computational nodes that need to be included in the biharmonic stencil correspond to those where $\psi<0$. By solving (11), this approach provides a continuous expression for $F$ over the entire domain while maintaining $V_{n}=F$ on the interface. We refer the reader to [58] for more details.

We now provide a simple test on our numerical scheme. As we show in Appendix A, the far-field boundary condition (8d) results in the rate of change of volume of a bubble evolving according to (8a)-(8d) to be $-4\pi$, independently of $\sigma$. Once a change in topology has occurred, the total rate of change of volume of all bubbles will remain $-4\pi$. To confirm this is consistent with our numerical solution, we consider the contraction rate of both a symmetric and asymmetric bubble, shown in Fig. 3, both of which undergo pinch-off. Fig. 4 shows the rate of decrease of volume $\dot{V}$, of both satellite bubbles, the exact rate of decrease of volume, and the rate of decrease of the satellite bubbles which form post pinch-off. For the symmetric bubble, Fig. 4$(a)$ indicates that the total rate of change of volume of the numerical solution compares well with the expected rate. Furthermore, the rate of decrease of the satellite bubbles are both $\dot{V}=-2\pi$, as expected. For the asymmetric bubble, Fig. 4$(b)$ shows that the total volume decreases at the expected (constant) rate both before and after pinch-off occurs.

We note that our model assumes the pressure of each bubble remains equalised between the two disconnected interfaces. In turn, this can lead to counter-intuitive volume changes. For example, if the inviscid bubble is neither injected or withdrawn ($Q=0$) such that $\dot{V}=0$, and if the bubble pinches-off into two satellite bubbles of unequal volume, as shown in Fig 3$(b)$ for example, then it is expected that the smaller bubble will contract to a point. In turn, this would result in the larger bubble expanding. However, as our primary motivation is to investigate the behaviour of solutions to (5a)-(5d) both approaching and recoiling from a pinch-off singularity, we expect that our model is still reasonable.

III.2 Pinch-off for a contracting bubble

III.2.1 Similarity exponents

We wish to study solutions like those presented in Fig. 3 in the neighbourhood of the point at which pinch-off occurs, $z=z_{0}$. If we set the time at which pinch-off occurs to be $t=t_{0}$, then we are particularly interested in the solution behaviour in the limit $t\rightarrow t_{0}^{-}$. It proves useful to use cylindrical polar coordinates so that $\phi=\phi(\rho,z,t)$ and the interface $\partial\Omega(t)$ is denoted by $\rho=h(z,t)$. In these coordinates, (5a)-(5c) become

[TABLE]

We seek similarity solutions to (12a)-(12c) of the form

[TABLE]

[TABLE]

In these similarity variables, the moving boundary $\rho=h(z,t)$ is described by $\xi=f(\zeta)$. By balancing the time-dependence on each side of (12a)-(12c), these equations become

[TABLE]

provided

[TABLE]

To formulate far-field conditions for $f(\zeta)$, we require that $\partial h/\partial t$ is finite away from $z=z_{0}$, which implies that $-f/3+\zeta f^{\prime}/3\rightarrow 0$ as $\zeta\rightarrow\pm\infty$ or, alternatively,

[TABLE]

By rescaling $\Phi$ by a further $\sigma^{2/3}$ and $f$, $\xi$, and $\zeta$ by $\sigma^{1/3}$, the surface tension parameter $\sigma$ can be completely eliminated from (15a)-(15c), suggesting that the pinch-off behaviour is universal and independent of all parameters. In terms of our dimensional quantities, the resulting similarity transformations would become

[TABLE]

We do not attempt to solve (15) and (17), but instead proceed to derive a more tractable system for the similarity solution below. For now it is worth emphasising that we are able to derive the similarity exponents (16) without any further analysis, which implies that (13) is a similarity solution of the first kind. Further, it is interesting to note the scaling exponent $\alpha=1/3$ was observed experimentally in [47] when 3He crystals pinched off; we discuss this issue further in Sec. V.

III.2.2 Similarity solution via an integral equation

In the following our approach is similar to that performed in Eggers and coworkers [15, 17] who study the break up of a bubble immersed in a fluid of low viscosity. In cylindrical polar coordinates, a solution to Laplace’s equation (12a) can be represented as the integral

[TABLE]

where the bubble has a length of $2L(t)$ and $C(z,t)$ is an unknown function to be determined. The kinematic boundary condition (12b) becomes

[TABLE]

while the dynamic boundary condition (12c) becomes

[TABLE]

Equations (20)-(21), which form a closed form system for $h(z,t)$ and $C(z,t)$, are an analogue of the system derived in [15, 17] for a bubble in a inviscid fluid.

As in Sec. III.2.1, we write $h$ out as a similarity solution as in (13)-(14), where again $t=t_{0}$ is the pinch-off time. We also write out $C(z,t)=(t_{0}-t)^{\gamma}D(\zeta)$. By balancing exponents, we find that provided $\alpha$, $\beta$ and $\gamma$ are given by (16), then (20)-(21) become

[TABLE]

[TABLE]

Note that for this similarity solution we have taken the limit $L/(t_{0}-t)^{1/3}\rightarrow\infty$. We solve (22)-(23) numerically for $f$ and $D$, as described below.

Before comparing the similarity solution from (22)-(23) with our numerical solutions to the full moving boundary problem, we proceed to make an approximation. We assume the main contribution to the integral in (22) is local around $\hat{\zeta}=\zeta$ and expand the integrand about this point. To leading order we find

[TABLE]

which suggests $D\approx f(f-\zeta f^{\prime})/6$. As a result, we have from (23) a single equation for $f(\zeta)$:

[TABLE]

In contrast to the problem of bubble pinch-off in an inviscid fluid, the lengths $\rho=h$ and $z$ in our Darcy flow problem scale the same way (both with the exponent $1/3$), so this type of analysis is not exact in the limit $t\rightarrow t_{0}^{-}$. However, as we now show, solutions to (25) with (17) still provide a very good approximation to the pinch-off profile, which is interesting.

We compare the numerical solution of (22)-(23) with the numerical solution to (25), shown in Fig. 5. Solutions are computed on the truncated domain $-40\leq\xi\leq 40$, which is discretised into equally spaced grid points with $\Delta\xi=0.1$. The derivatives and integrals are approximated using central differences and the trapezoidal rule, respectively. This approach leads to a nonlinear system of algebraic equations which is solved using the Newton-Raphson method. We find that the gradient of $f$ from the solution to (22)-(23) is slightly steeper than the solution to (25) in the far field. However, around the pinch-region where $\xi=0$, the two solutions are indistinguishable (at this scale). This comparison suggests that the solution to (25) works effectively as well as the solution to (22) and (23) for describing the behaviour of bubbles approaching a pinch-off event. This observation is interesting since the spatial rescaling in the $\rho$ and $z$-directions is the same, suggesting that the neck is not slender except very close to point at which $\partial h/\partial z=0$.

The farfield boundary condition (17) implies the linear behaviour $f\sim\pm A\xi$ as $\xi\to\pm\infty$ for some unknown $A$ such that the free surface takes the shape of a double cone at the point of pinch-off. We approximate $A$ from the numerical solution to (22) and (23) to be $A\approx 0.97$ which gives a cone angle of $\theta\approx 44^{\circ}$. In comparison, the numerical solution to (25) gives $A\approx 0.94$ resulting in a slightly smaller angle of $\theta=43^{\circ}$. We note that this cone angle is larger than the value determined experimentally by Ishiguro et al. [47], who approximated this angle from a video frame of the experiment just after pinch-off had occurred and found $\theta\approx 35^{\circ}$.

III.2.3 Numerical solution to full moving boundary problem

We now present numerical solutions to the full moving boundary problem (5a)-(5d), produced by our numerical scheme described in Sec. III.1. We begin by checking whether the similarity exponent of $\alpha=1/3$ is consistent with the full solution. To do so, we solve (8a)-(8d) for the case in which the interface is initially an asymmetric dumbbell (see Fig. 3) for different values of the surface tension parameter $\sigma$. While increasing $\sigma$ decreases the time $t_{0}$ at which pinch-off occurs, the similarity exponent $\alpha$ is independent of the parameters of the model. Fig. 6$(a)$ compares the minimum neck radius of each of these solutions solution as a function of the time to pinch-off, $t_{0}-t$, where $\sigma=1$, (blue), 0.5 (red), and 0.1 (yellow), while Fig. 6$(b)$ shows the corresponding $\log$-$\log$ plot. By taking a line of best fit (black dashed lines) for each value of $\sigma$, we compute $\alpha=0.33$ to two decimal places, which is in good agreement with the theoretical result $\alpha=1/3$ derived in Sec. III.2.1.

We also compare the interfacial profiles of the numerical solution to (8a)-(8d) with the similarity solution (22)-(23). Fig. 7 shows the evolution of the solution to (8a)-(8d) for both a symmetric and asymmetric initial conditions (see Fig. 3) up to $t=t_{0}$. The inserts compare the numerical solution of (8a)-(8d) scaled according to (13)-(14) with the numerical solution to (22)-(23). For both the symmetric and asymmetric initial conditions, around the pinch region we observe good agreement between the similarity solution and the full numerical solution. Away from the pinch region where $|\xi|$ is large, there is noticeable deviation between the numerical solution to (8a)-(8d) and the solution to (22)-(23), at least for the asymmetric case. As (22)-(23) was derived in the limit that pinch-off is approached, these deviations are to be expected.

III.3 Recoil

In Sec. III.2, we investigated the behaviour of the bubble whose boundary’s velocity blows up as it approaches pinch-off due to the unbounded increase in interfacial curvature. Once pinch-off has occurred, the two newly formed satellite bubbles will rapidly recoil due to the large local curvature at each of the tips. While many studies concerned with the pinch-off of viscous or inviscid fluids have focused on the behaviour of the interface approaching pinch-off, fewer studies consider the behaviour of post break-up [59, 60, 61]. As a starting point, we expect the argument in Sec. III.2.1 to still hold post break-up, leading to the same similarity exponents, namely $\alpha=1/3$ and $\beta=1/3$, as was computed from the pre pinch-off solution.

To confirm that our numerical solution produces the scaling exponent of $\beta=1/3$, we consider the evolution of an asymmetric dumbbell bubble (see Fig. 3) for $t>t_{0}$. After pinch-off, we compute the distance between the point where pinch-off occurs and the tip of the two retracting satellite bubbles as a function of the time. Fig. 8$(a)$ shows the evolution of this distance for both the left and right satellite bubbles, while Fig. 8$(b)$ is the corresponding plot on a $\log$-$\log$ scale. By taking a line of best fit, we find an approximation of the similarity exponent to be $\beta=0.33$ to two decimal plates for both bubbles, which is consistent with the theoretical prediction. We repeat this procedure for a range of values of $\sigma$ (not shown) and, as was found in Sec. III.2.3, our numerical results are in good agreement with the theoretical prediction of $\beta=1/3$.

Just after pinch off, we expect the interface to also act in a self-similar fashion. In this case, however, the representation (19) is unlikely to be valid [17], and the approximation (24) no longer true, as $f^{\prime}(\xi)\to\infty$ where the interface intersects the $\xi$-axis. Due to this difficulty we do not solve the similarity problem here. While we do not provide a numerical solution of the asymptotic equations post break-up, we perform a numerical test to confirm that solutions to (8a)-(8d) are self-similar in the limit. We redefine our self-similar solution after pinch-off to be

[TABLE]

noting that now $t>t_{0}$. In Fig. 9, we plot the numerical solution to (8a)-(8d) scaled according to (26) for times where $t-t_{0}\ll 1$ for the $(a)$ symmetric and $(b)$ asymmetric dumbbell initial conditions. For both cases, Fig. 9 suggests that the interfaces do indeed converge to a self-similar profiles around the pinch region.

IV Translating bubbles in a cylindrical tube

We now turn our attention to a situation in which an inviscid bubble is translating from left to right in a cylindrical tube of unit diameter, as illustrated in Fig. 1$(b)$-$(c)$. In Fig. 1$(b)$, the viscous fluid fills the domain to the right of the bubble interface, while in Fig. 1$(c)$ the bubble is finite and the surrounding viscous fluid extends infinitely far in both directions. In cylindrical polar coordinates with the $z$-axis running down the centre of the cylinder, velocity potential is $\phi=\phi(\rho,z,t)$ and the interface is denoted by $\rho=h(z,t)$. We impose a no-flux boundary condition on the cylinder wall so that $\partial\phi/\partial\rho=0$ on $\rho=1/2$. Further, our flow is driven by a flux of fluid in the far field so that $\partial\phi/\partial z\sim 1$ as $z\to\infty$. It is straightforward to adapt the numerical scheme presented in Sec. III.1 to this new geometry.

IV.1 Development of a Saffman-Taylor-like finger

We consider here the geometry in Fig. 1$(b)$ where the inviscid fluid is displacing the viscous fluid. The Saffman-Taylor instability causes a small perturbation of the interface to grow and eventually a single long finger of inviscid fluid will develop and propagate into the viscous fluid. To illustrate this progression, we show in Fig. 10 a solution for which the initial condition is a flat interface with a small sinusoidal perturbation. The evolution to a long finger is clear.

This example is an axisymmetric analogue of the well-studied Hele-Shaw problem for which a Saffman-Taylor finger develops in a rectangular channel [25]. As is well known in that two-dimensional problem, both experimentally and in time-dependent numerical simulations, a single finger develops whose width $\lambda$ appears to be slightly more than half of the channel diameter provided the capillary number is sufficiently high [62, 63]. Assuming a travelling wave solution with a finger moving at speed $U$, numerical studies show there is a countably infinite number of possible finger widths $\lambda$ for each fixed value of the surface tension $\sigma$ [24, 64, 65]. The lowest of these finger widths has been shown to correspond to a stable solution while the remaining solutions are unstable [66]. Each family of solutions has the asymptotic behaviours $\lambda\rightarrow 1/2^{+}$ as $\sigma\rightarrow 0$ and $\lambda\rightarrow 1^{-}$ as $\sigma\rightarrow\infty$ [67, 68, 69, 70].

Returning to the axisymmetric version of this problem with which we are concerned with here, there are numerical results for the travelling wave problem by Levine & Tu [13], who show that much of the behaviour is qualitatively similar to the Hele-Shaw case, with multiple solutions for fixed values of the surface tension $\sigma$. However, a key difference is that the solution branches do not all tend towards $\lambda=1/2$ as $\sigma$ decreases, but instead solutions only exist down to a minimum value of $\sigma$ at which point pairs of branches merge together (this type of merging of pairs of branches also occurs for the Hele-Shaw wedge problem studied by Ben Amar & Combescot [71, 72, 73]).

Our contribution here is to study the full time-dependent axisymmetric problem which, as far as we know, has not been documented before. To gather data equivalent to that provided in [13], we must compute the transient solution for a sufficient amount of time such that the finger which develops evolves at a constant shape and speed. Of course, this approach only allows us to track the stable branch, as the solution will not evolve to an unstable travelling wave solution.

We show in Fig. 11 the results of our calculations with various values of the finger width $\lambda$ plotted against surface tension $\sigma$. The smallest value of finger width we observe is $\lambda\approx 0.6$, which is similar to the minimum value observed by Levine & Tu [13]. Our time-dependent calculations for smaller values of surface tension do not lead to single steadily propagating fingers, but instead a type of axially symmetric tip-splitting occurs, where a dimple forms at the front of the bubble and subsequently the bubble takes the shape of an off-axis finger rotated around the $z$-axis. We speculate that this bubble shape would not be observed experimentally as it is likely unstable to non-axisymmetric perturbations.

IV.2 Pinch-off of a Saffman-Taylor-like finger

The simulations discussed in Sec. IV.1 were initiated with a small perturbation of flat interface. For very different non-convex initial conditions which involve a thin neck, the bubble can undergo a pinch-off event which is similar to that described in Sec. III.2. One such example is provided in Fig. 12, illustrating how the very large azimuthal curvature in the neck region ‘pulls’ the interface towards pinch-off even though the volume is increasing. These singularities do not develop in the classical two-dimensional Hele-Shaw problem as there is no increasingly large component of curvature in that case. As such, this behaviour is another example of a difference between the axially symmetric problem we consider here in this paper and the Hele-Shaw version.

We now demonstrate that the pinch-off singularity that occurs in this cylindrical geometry is the same as that presented in Sec. III.2 for a shrinking bubble. Following the same methodology, we estimate $\alpha$ from the full numerical solution by taking a line of best fit of the minimum neck radius $h_{\mathrm{min}}$ as a function of the time to pinch-off, $t_{0}-t$, on a $\log$-$\log$ scale, as shown in Fig. 13$(a)$. We find a value of $\alpha=0.33$ to two decimal places, which is again consistent with both the similarity problem derived in Sec. III.2.1 and the numerical solution of the contracting bubble presented in Sec. III.2.3. Furthermore, the insert compares the full numerical solution scaled according to (13) with the solution to (22)-(23). As we observed in Fig. 7, around the pinch region where $|\xi|$ is small we find good agreement between both the similarity solution and the full numerical solution. Regarding the post pinch-off behaviour, Fig. 13$(b)$ shows a $\log$-$\log$ plot of the distance between point where pinch-off occurs and the tip of the recoiling bubble as a function of time after pinch-off, $t-t_{0}$. Again, we estimate a similarity exponent consistent with the theoretical prediction.

IV.3 Evolution of a Taylor-Saffman-like bubble

In Sec. IV.2 we discussed how an interface may break up to produce a bubble of finite volume that transverses through a cylindrical tube. When such a bubble propagates along the channel and eventually settles into a steady motion with constant shape and speed $U$, it can be thought of as an axially symmetric version of a Taylor-Saffman bubble [26]. By choosing the initial condition of the bubble to be a prolate spheroid with varying aspect ratios, we are able to plot the steady-state speed $U$ versus surface tension $\sigma$ for a fixed volume of these Taylor-Saffman-like bubbles, as in Fig. 14. These plots show the same qualitative behaviour as the well-studied two-dimensional case [27, 74]. Namely, for small bubbles ($\blacktriangle$ and $\blacktriangledown$), the steady state shape of the bubble will tend to a sphere as $\sigma$ increases, whose speed is dependent upon the radius of this sphere. For larger bubbles whose size is too large to allow a spherical shape in the channel ($\blacktriangleleft$ and $\blacktriangleright$), their speed appears to be approaching unity as $\sigma$ increases, although it is not clear whether there is a maximum value of $\sigma$ for which solutions on this branch exist. A difference between our axially symmetric bubbles and the two-dimensional Hele-Shaw problem is that in our case, each solution branch approaches $U=3^{-}$ in the limit $\sigma\rightarrow 0$, while in the two-dimensional case, the branches select $U=2^{-}$ in this limit.

As just mentioned, while small steady bubbles will tend to spheres for large surface tension, the limiting shape of large steady bubbles as surface tension is increased is less obvious. We compare the numerical solution computed using different values of the surface tension parameter (at a sufficiently large time such that the bubble’s speed and shape are essentially constant) in Fig. 15. For low surface tension (Fig. 15$(a)$), we find that the bubble tends to a long thin shape. As surface tension increases (Fig. 15$(b)$), we find that the sides of the bubble are no longer parallel but instead the bubble’s trailing end is wider. For large values of surface tension (Fig. 15$(c)$ and $(d)$), we find the length of the bubble decreases and the width of the bubble is nearly as large as the radius of the cylinder.

The shape of the bubble in Fig. 15$(b)$ is similar to those of Taylor bubbles in cylindrical tubes. This problem, where a large bubble is moving steadily in a tube otherwise filled with viscous fluid (in contrast to our problem, where the tube is filled with a porous medium saturated with viscous fluid), has an extensive history going back to Bretherton [75] and Taylor [76] and is still attracting plenty of interest (see [77, 78], for example). Our brief study in this section suggests that it is worth investigating further the links between our Taylor-Saffman-like bubbles and the well-studied Taylor bubbles.

V Discussion

In this paper, we conducted a numerical and theoretical study investigating the behaviour of solutions to a one-phase Darcy model describing the evolution of a inviscid bubble within a fully saturated homogeneous porous medium. The model is a moving boundary problem which is an analogue of that for the one-phase Hele-Shaw problem, including a dynamic condition (3) with pressure across the interface related to the curvature $\kappa$ via the surface tension. This type of model is analysed and discussed in [7, 8, 9, 6, 10, 11, 12, 13, 3, 4, 5, 14], for example. An important difference between our axially symmetric model and the Hele-Shaw version is that $\kappa$ represents the mean curvature of the surface while in the Hele-Shaw model it represents curvature in the plane. This difference drives new behaviour not witnessed in the Hele-Shaw problem, which is worth recording.

A key focus of our paper is to study how a bubble undergoes a pinch-off singularity. In Sec. III.1, we summarised a flexible numerical scheme based on the level set method that is able to accurately solve (5a)-(5d) both before and after a pinch-off singularity has developed. Using this scheme, we have shown that the self-similar solution derived by applying an integral equation (motivated by [15]) matches well with full numerical solutions to (5a)-(5d). This was done for both contracting bubbles in radial geometry (Sec. III) and translating bubbles in cylindrical channel geometry (Sec. IV).

It is worth mentioning the problem of an expanding bubble due to injecting an inviscid fluid at a point, as illustrated by Fig. 1$(d)$. This is a very well studied geometry in the two-dimensional Hele-Shaw case, as it acts as a prototype problem for viscous fingering and interfacial pattern formation [79, 80, 81]. In three dimensions, the interfacial dynamics are more complicated due to the extra degree of freedom [9, 12] and the possibility of fingers pinching off the main bubble in the same manner as we have been studying in this paper. We illustrate this point via an example shown in Fig. 16. Here we solve (5a)-(5d) with the positive sign in (5d) in order to simulate a bubble expanding. For the particular initial condition used in Fig. 16, we see a finger propagate along the positive $z$-axis which eventually breaks off the main bubble. While we have not included the details here, the accompanying singularity in curvature is of the same form, both before and after pinch-off, as that considered in Sections III and IV. Of course, this example is partly contrived as it involves an axially symmetric geometry that does not illustrate the full three-dimensional nature of viscous fingering. It does, however, demonstrate how singularities in curvature develop for expanding bubbles, a feature which is not present in two dimensions.

While our main physical motivation has been to study the evolution of a bubble in a porous medium, as noted in Sec. II.2, (5a)-(5d) can also be interpreted as describing the melting/freezing of a crystal dendrite. The pinch-off of crystals has been experimentally studied in [47], where it was shown that the neck of the crystal both approaches and recoils from pinch-off with the scaling exponent $\alpha=1/3$. To make analytical progress, the authors argued that the solid-melt interface can be approximated via mean curvature flow, $v_{n}=-\kappa$, which was shown to have a similarity exponent of $\alpha=1/2$. The authors offer the explanation for this discrepancy that the mean curvature flow is slow to converge to self-similar behaviour, and the experimental results were not accurate enough to reflect this. It is interesting that our model (5a)-(5d), which is the large Stefan number version of (6) with (5b)-(5d), gives the same similarity scaling.

Apart from studying pinch-off singularities, our other key focus has been on the time-dependent development of axially symmetric analogues of Saffman-Taylor fingers and Taylor-Saffman bubbles in a cylindrical channel, highlighting the differences and similarities between the axially symmetric and two-dimensional problems. One noteworthy property of the finger problem is that families of steadily propagating axially symmetric fingers merge for finite values of surface tension in a similar manner to that observed in the Hele-Shaw wedge problem [71]. Our numerical results indicate that the speed of bubbles (of any volume) will tend towards $3^{-}$ in the limit that $\sigma\to 0$, which differs from Taylor-Saffman bubbles in a two-dimensional Hele-Shaw cell whose speed is known to tend towards $2^{-}$ [82, 83, 27, 74, 84]. More accurate calculations of the data in Fig. 14 and shapes of steady-state axially symmetric Taylor-Saffman-like bubbles could be obtained by moving to the reference frame of the steadily moving bubbles and applying a boundary integral method. The accompanying selection problem for $U=3$ would provide an interesting challenge requiring techniques in exponential asymptotics, especially given the lack of exact solution for the leading order problem.

Acknowledgements

SWM and LCM acknowledge the support of the Australian Research Council Discovery Project DP140100933. We are grateful to the anonymous referees for their helpful feedback.

Appendix A Prescribed rate of volume decrease

We now show that the form of the far-field boundary condition (5d) results in the rate of change of volume of a bubble evolving according to (5a)-(5d) to be $\pm 4\pi$. The rate of change of volume of a hypersurface whose boundary moves with normal velocity $v_{n}$ is

[TABLE]

which, using the kinematic boundary condition (5b), gives

[TABLE]

Supposing that the viscous fluid occupies the region \mbox{\boldmathx}\in B\backslash\Omega(t), where $B$ is a sphere of radius $R$, then

[TABLE]

which comes about via the divergence theorem. Setting $\nabla^{2}\phi=0$ and substituting (28) into (29) gives

[TABLE]

We note that if a bubble were to undergo a change in topology, the rate of change of the total volume will remain $\pm 4\pi$. However, each individual bubble may not expand or contract at a constant rate. We further investigate the evolution of bubbles after pinch-off has occurred in Sec. III.1.