2. Film bursting at an air–liquid interface: experimental method

We study the three-phase configuration experimentally by inflating an oil drop (‘s’ for ‘surroundings’ in figure 1c) at a water–air free interface and capturing the retraction of the water film ( f in figure 1c). The schematic of the experimental set-up is shown in figure 2(a). A plastic box (from the Dutch store ‘Bodemschat’) of dimensions $25\ {\rm mm} \times 25\ {\rm mm} \times 15\ {\rm mm}$ (length $\times$ width $\times$ height) filled with purified water (Milli-Q) was used as the liquid bath for most of the experiments. To study the effect of the viscosity of the retracting film, the water in the bath was replaced by glycerol–water (glycerol from Sigma-Aldrich) mixtures (concentrations in the range 50 %–70 % by weight) for some experiments. A dispensing needle (${\rm inner\ diameter} = 0.41\ {\rm mm}$, HSW Fine-Ject) was submerged within the bath such that its dispensing end was at a depth of $2.4\ {\rm mm}$ from the free surface (depth kept constant during all experiments). A silicone oil (Wacker) droplet was created at the tip of the needle by connecting it to an oil-filled plastic syringe ($5\,\,{\rm ml}$, Braun Injekt) via a flexible plastic polyetheretherketone (PEEK) tubing (Upchurch Scientific). The oil flow rate was maintained at $0.05\,\,{\rm ml}\,{\rm min}^{-1}$ with the help of a syringe infusion pump (Harvard Apparatus). In the present experiments, silicone oils of different viscosities were used, and their densities ($\rho _{s}$), kinematic ($\nu _{s}$) and dynamic viscosities ($\eta _{s}$) are listed in table 1. It is to be noted that for only the AK 0.65 oil ($\eta _{s} = 4.94 \times 10^{-4}\,{\rm Pa}\,{\rm s}$), the drop is less viscous than the water film ($\eta _{f} = 8.9 \times 10^{-4}\,{\rm Pa}\,{\rm s}$), while for all the other oils, the film is less viscous. The oil–water interfacial tension ($\gamma _{sf}$) was considered to be $0.040\,{\rm N}\,{\rm m}^{-1}$ (Peters & Arabali Reference Peters and Arabali2013).

Figure 2. (a) Schematic of the experimental set-up. (b) Typical time-lapsed experimental snapshots of the film rupture and the subsequent retraction process ($\nu _{s} = 10\,{\rm cSt}$). The time instant $t$ = $t_{0}$ denotes the first frame where rupture (indicated by the white arrow) is discernible.

Table 1. Salient properties of the silicone oils used in the present work.

The drop volume was increased by a slow infusion using the syringe pump. The needle depth below the free surface was chosen such that the drop remained anchored to the needle during inflation. As a result, the water film right above the oil droplet progressively thins with increasing volume of the drop. Below a certain thickness, the film ruptures due to van der Waals forces (Vaynblat, Lister & Witelski Reference Vaynblat, Lister and Witelski2001), and subsequently retracts into the bath. This situation is analogous to the rupture and retraction of a water film sandwiched between air and a viscous oil droplet. This is also a configuration that is flipped vertically as compared with the early-time scenario studied by Feng et al. (Reference Feng, Muradoglu, Kim, Ault and Stone2016). Another key difference is that we ensure negligible vertical velocity at the point of rupture of the water film, making this scenario ideal for studying three-phase Taylor–Culick retractions.

High-speed imaging of this rupture and retraction phenomena was performed at 50 000 f.p.s. (frames per second) for the lower viscosity oils and at 10 000 f.p.s. for the higher viscosity ones, with a $2.5\,\mathrm {\mu }{\rm s}^{-1}$ exposure time, by a high-speed camera (Fastcam Nova S12, Photron) connected to a macro lens (DG Macro $105\ {\rm mm}$, Sigma) with $64\ {\rm mm}$ of lens extender (Kenko). The camera was pointed at a plane mirror (Thorlabs) inclined at $45^{\circ }$ to the horizontal to capture the top view of the retraction phenomenon (figure 2a), while the experiments were illuminated from the top by a light-emitting diode light source (KL 2500 LED, Schott). A typical bursting event is shown in figure 2(b), where time $t=t_{0}$ indicates the instant when rupture is optically discernible. With increasing time, the size of the hole formed due to rupture increases as the film retracts. The phenomenological observations shown in figure 2(b) will be discussed in detail in § 3. The captured images were then further analysed using the open-source software FIJI (Schindelin, Aganda-Carreras & Frise Reference Schindelin, Aganda-Carreras and Frise2012) and an in-house OpenCV-based Python script to obtain quantitative information presented in the following sections.

3. Film bursting at an air–liquid interface: experimental results

The rupture and retraction of a water film on the surface of an oil drop of $\eta _{s}=4.94\times ~10^{-4}\,{\rm Pa}\,{\rm s}$ is shown in figure 3(a) (and supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.671). The timestamps indicate ($t - t_{0}$), where $t_{0}$ is the time instant when rupture is optically discernible, and $t$ is the current time. As mentioned earlier, in this particular case, the water film is more viscous than the oil. It is to be noted that in the present experiments, we could not precisely control the location of rupture as it was sensitive to experimental noise (see § 4.2 of Villermaux (Reference Villermaux2020)). Hence, the rupture in the present experiments did not always occur at the apex of the thinning film. Such behaviour was also observed in other similar experiments of film rupture (Oldenziel, Delfos & Westerweel Reference Oldenziel, Delfos and Westerweel2012; de Maleprade, Clanet & Quéré Reference de Maleprade, Clanet and Quéré2016). The rupture location may also be determined by a ‘prehole’ formation, also observed by Vernay et al. (Reference Vernay, Ramos and Ligoure2015) for the bursting of emulsion-based liquid sheets. In their work, Vernay et al. (Reference Vernay, Ramos and Ligoure2015) show that the presence of emulsion oil droplets at the air–water interface results in lowering of the local interfacial tension, leading to Marangoni flows away from that location. This flow leads to a local thinning of the film, which ultimately ruptures at that location. They also report that the prehole formation always precedes rupture in their experiments. In the present experiments, the water surface is never pristine and always contains small impurities (which are practically unavoidable). It is possible that these impurities might have reduced the local surface tension, resulting in a similar Marangoni flow leading to a prehole. For the discussion on the origin of the hole nucleation, we also refer to Lohse & Villermaux (Reference Lohse and Villermaux2020). In any case, upon rupture, a circular hole is formed in the film, which grows radially in time. Therefore, the oil bounded by the periphery of the hole gets into contact with air and not with the water film. It is also noticeable that the edge of the retracting film forms a thick rim – an observation also made for the retraction of liquid films in air (Pandit & Davidson Reference Pandit and Davidson1990; Brenner & Gueyffier Reference Brenner and Gueyffier1999; Sünderhauf et al. Reference Sünderhauf, Raszillier and Durst2002). The presence of the rim can be qualitatively surmised from the experimental snapshots (figures 2b and 3a,b), where the change in the curvature of the film downstream of the rim introduces a difference in the colour intensity. As the hole increases in size (or as the film retracts further), this rim also becomes thicker. Finally, since silicone oil prefers to spread on water (Li et al. Reference Li, Diddens, Segers, Wijshoff, Versluis and Lohse2020), the retraction process ceases when the oil droplet has completely spread on water, thus creating a macroscopic film whose thickness is controlled by volume conservation and thermodynamics (de Gennes, Brochard-Wyart & Quéré Reference de Gennes, Brochard-Wyart and Quéré2004).

Figure 3. Time-lapsed snapshots of the postrupture retraction of water films for different oil viscosities: (a) $\nu _{s}$ = 0.65 cSt; (b) $\nu _{s}$ = 10 cSt; and (c) $\nu _{s}$ = 100 cSt. The rupture location is denoted by the white arrows, while the time stamps indicate the time since rupture is first observed (i.e. $t - t_{0}$). Also see supplementary movie 1.

When the viscosity of the oil phase is increased to $\eta _{s} = 9.30 \times 10^{-3}\,{\rm Pa}\,{\rm s}$, the hole opening (or film retraction) dynamics (as seen in figure 3b and supplementary movie 1) are qualitatively similar to that for the lower viscosity described before (figure 3a). Here also, the film forms a thick rim at its retracting edge. Nonetheless, an increase in the oil's viscosity decreases the retraction speed, as indicated by the timestamps (corresponding to $t - t_{0}$). This behaviour is expected since the physical situation is analogous to a retracting water film shearing the free surface of viscous oil: increasing $\eta _{s}$ increases the resistance to shearing, which in turn makes the retraction process slower.

Furthermore, the experimental snapshots show that the oil–air–water contact line exhibits corrugations during the retraction process, and fine streams of droplets are released from these corrugations. Similar observations were also made for film retraction in the two-phase configuration (Reyssat & Quéré Reference Reyssat and Quéré2006; Oldenziel et al. Reference Oldenziel, Delfos and Westerweel2012) and during the rupture of the intermediate film when a drop coalesces with a pool of the same liquid in the presence of an external medium (Aryafar & Kavehpour Reference Aryafar and Kavehpour2008; Kavehpour Reference Kavehpour2015). The nature of the corrugations is reminiscent of the sharp tips observed during selective withdrawal (Cohen & Nagel Reference Cohen and Nagel2002; Courrech du Pont & Eggers Reference Courrech du Pont and Eggers2006, Reference Courrech du Pont and Eggers2020) or tip streaming (Montanero & Gañán Calvo Reference Montanero and Gañán Calvo2020). In a frame of reference comoving with the rim, the film sees highly viscous oil being aspirated away from it, resulting in the formation of the sharp tips. Indeed, such a mechanism was also hinted at by Reyssat & Quéré (Reference Reyssat and Quéré2006) for the instabilities observed in their experiments for film retraction in the two-phase configuration. Tseng & Prosperetti (Reference Tseng and Prosperetti2015) showed that such instabilities are formed due to the local convergence of streamlines in the neighbourhood of a zero-vorticity point or line on the interface. However, a detailed and quantitative investigation of the formation and subsequent breakup of these liquid tips is beyond the scope of the present work.

For an even higher viscosity of the oil phase ($\eta _{s} = 9.60 \times 10^{-2}\,{\rm Pa}\,{\rm s}$; see figure 3c and supplementary movie 1), the retraction of the ruptured water film is further slowed down (as evident from the timestamps in figure 3c). Furthermore, the retracting edge also does not possess a thick rim. This observation is similar to the case of Brenner & Gueyffier (Reference Brenner and Gueyffier1999) for the retraction of viscous films in air, where films of higher viscosity do not form a rim. Moreover, although the expanding holes for the lower $\eta _{s}$ cases (as shown in figure 3a,b) are almost circular, the one for the high viscosity case shown in figure 3(c) is highly asymmetric. This asymmetry can be attributed to the location of the rupture not being at the film's apex. Since the rupture is happening at an off-apex location, the film thickness at the location of rupture is not spatially uniform due to the curvature of the oil droplet. Hence, the retraction velocity is faster on the part of the film which has a lower thickness. Presumably, this effect is more pronounced when the overall film retraction dynamics are slower, as is the case for the experiments shown in figure 3(c). To confirm this hypothesis, one requires high-resolution measurements of the spatial variation of the film thickness, which is challenging in the present experiments (further discussed in § 7). The corrugations at the oil–air–water contact line are also observed in this case. However, since the retraction velocity itself is considerably smaller than for the case shown in figure 3(b) (see figure 4b for specific values), the tips are not as sharp, and no droplet streams are observed. To quantify the retraction dynamics, we measure the hole opening radius from each snapshot captured using the high-speed camera. For each experimental snapshot, the area of the hole $A(t)$ is measured, and subsequently an equivalent hole opening radius $R(t)$ is calculated as $A(t) = {\rm \pi}(R(t))^{2}$. A typical measurement from the optical images is depicted in the inset of figure 4(a). The temporal variation of the measured hole radius, $R$, is shown in figure 4(a). The time instant corresponding to the first frame in which rupture is optically discernible is denoted by $t_{0}$. Each datapoint in figure 4(a) denotes the mean of measurements from five independent experiments, and the error bars correspond to $\pm$ one standard deviation. In the present work, we focus on the early moments following rupture, as indicated by the red rectangle in figure 4(a). Zooming into this early-time regime, as shown in figure 4(b), it is observed that $R(t)$ varies linearly with time (as evident from the lines denoting linear fits in figure 4b). This variation indicates that the retraction velocity $v_f$ (= ${\rm d}R/{\rm d}t$), given by the slopes of the linear fits, is constant for each viscosity. This is reminiscent of the constant rupture velocity also observed for the classical (figure 1a) and two-phase (figure 1b) Taylor–Culick configurations. Furthermore, it is also observed that with increasing $\eta _{s}$ (or $\nu _{s}$), the slope of the linear fits (hence $v_f$) decreases, as expected from the qualitative observations reported in figure 3.

Figure 4. (a) Temporal evolution of the retraction radius ($R$) for oils of different kinematic viscosity ($\nu _{s}$); a typical measurement is shown in the snapshot in the inset. (b) At early times, red rectangle in panel (a), $R$ varies linearly with time; the discrete datapoints are experimental measurements and the lines are linear fits. (c) Variation of dewetting velocity ($v_f$) with the dynamic viscosity of the oil phase ($\eta _{s}$); the discrete datapoints are experimental measurements and the line represents $v_f \sim \eta _{s}^{-1/2}$.

The variation of $v_f$ with $\eta _{s}$ is shown in figure 4(c). The typical retraction velocities are $O(1\,{\rm m}\,{\rm s}^{-1})$. A decreasing $v_f$ with increasing $\eta _{s}$ is observed. Furthermore, for the cases where the oil is more viscous than water ($\eta _{f} = 8.9 \times 10^{-4}\,{\rm Pa}\,{\rm s}$), the retraction velocity varies as

(3.1)\begin{equation} v_{f} \sim \frac{1}{\eta_{s}^{1/2}}, \end{equation}

as evident from the line in figure 4(c). This is a weaker dependence as compared with the expected $1/\eta _{s}$ variation observed for retraction in the two-phase configuration (Martin et al. Reference Martin, Buguin and Brochard-Wyart1994; Eri & Okumura Reference Eri and Okumura2010). We will attempt to explain the scalings for the two-phase and three-phase configurations in § 7. Furthermore, the reason for not fitting the datapoint for the case where the oil is less viscous than the water film in figure 4(c) will also be addressed therein.

4. Numerical framework

4.1. Governing equations

In this section, we discuss the governing equations that describe the retraction of a ruptured liquid film in the three configurations we study in this paper, namely, the classical, two-phase and three-phase Taylor–Culick retractions. We perform axisymmetric direct numerical simulations using the free-software VoF program, Basilisk C (Popinet & Collaborators Reference Popinet2013–2022; Sanjay Reference Sanjay2021b), which uses the one-fluid approximation (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011) to solve the continuity and Navier–Stokes equations,

(4.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v} = 0, \end{gather}$$

(4.2)$$\begin{gather}\rho\left(\frac{\partial\boldsymbol{v}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{vv}\right)\right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(2\eta\boldsymbol{\mathcal{D}}\right) + \boldsymbol{f}_\gamma, \end{gather}$$

where, $\boldsymbol {v}$ and $p$ are the velocity vector and pressure fields, respectively, $\eta$ the viscosity of the fluid and $t$ denotes time. Furthermore, $\boldsymbol {\mathcal {D}}$ is the symmetric part of the velocity gradient tensor $(\boldsymbol {\mathcal {D}} = (\boldsymbol {\nabla } \boldsymbol {v} + (\boldsymbol {\nabla } \boldsymbol {v})^{\text {T}})/2)$, and $\boldsymbol {f}_\gamma$ the singular surface tension force needed in the one-fluid approximation to comply with the dynamic boundary condition at the interfaces (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992).

4.2. Non-dimensionalization of the governing equations

We non-dimensionalize the governing equations by using the inertiocapillary velocity scale $v_\gamma$, the thickness of the film $h_0$ and the capillary pressure $p_\gamma$. These scales also define the characteristic inertiocapillary time, $\tau _\gamma$, as

(4.3a–c)\begin{equation} \tau_\gamma = \frac{h_0}{v_\gamma} = \sqrt{\frac{\rho_fh_0^3}{2\gamma_{sf}}}, \quad v_\gamma = \sqrt{\frac{2\gamma_{sf}}{\rho_f h_{0}}}, \quad p_\gamma = \frac{2\gamma_{sf}}{h_0}. \end{equation}

Here, $\gamma _{sf}$ is the surface tension coefficient between the film ( f) and the surrounding (s) medium, $\rho _f$ the film density and $h_0$ its thickness. The dimensionless form of the Navier–Stokes equation (4.2) is

(4.4)\begin{equation} \tilde{\rho}\left(\frac{\partial\tilde{\boldsymbol{v}}}{\partial \tilde{t}} + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left(\tilde{\boldsymbol{v}}\tilde{\boldsymbol{v}}\right)\right) ={-}\tilde{\boldsymbol{\nabla}} p + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left(2Oh\tilde{\boldsymbol{\mathcal{D}}}\right) + \tilde{\boldsymbol{f}}_\gamma, \end{equation}

where the expressions for the Ohnesorge number ($Oh$; the ratio of viscocapillary to inertiocapillary time scales), the dimensionless density ($\tilde {\rho }$) and the singular surface tension force ($\tilde {\boldsymbol {f}}$) depend on the specific configurations that we discuss below.

4.2.1. Two-phase Taylor–Culick configuration

In this configuration, a liquid film ( f) retracts in a viscous surrounding (s) medium (figure 5a). We use the VoF tracer $\varPsi$ to differentiate between the film ($\varPsi =1$) and the surroundings ($\varPsi = 0$), which follows the VoF scalar advection equation,

(4.5)\begin{equation} \left(\frac{\partial}{\partial\tilde{t}} + \tilde{\boldsymbol{v}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}\right)\varPsi = 0. \end{equation}

Furthermore, the singular surface tension force is given by (Brackbill et al. Reference Brackbill, Kothe and Zemach1992)

(4.6)\begin{equation} \tilde{\boldsymbol{f}}_\gamma \approx \left(\tilde{\kappa}/2\right)\tilde{\boldsymbol{\nabla}}\varPsi, \end{equation}

where the curvature $\kappa$ is calculated using the height-function approach (Popinet Reference Popinet2009). We follow the same sign convention as Tryggvason et al. (Reference Tryggvason, Scardovelli and Zaleski2011, page 33): the curvature is positive if the interface folds towards its normal $\boldsymbol {\hat {n}}$, i.e. $\kappa = - \boldsymbol {\nabla } \boldsymbol{\cdot } \boldsymbol {n}$. Note that the surface tension scheme in Basilisk C is explicit in time. So, we restrict the maximum time step as the characteristic inertiocapillary time based on the wavelength of the smallest capillary wave. Additionally, the density of the film is the same as that of the surroundings, giving $\tilde {\rho } = 1$. Lastly, the Ohnesorge number ($Oh$) is given by

(4.7)\begin{equation} Oh = \varPsi {Oh}_{f} + \left(1-\varPsi\right){Oh}_{s}, \end{equation}

where

(4.8a,b)\begin{equation} {Oh}_{f} = \frac{\eta_f}{\sqrt{\rho_f\left(2\gamma_{sf}\right)h_0}} \quad \text{and} \quad {Oh}_{s} = \frac{\eta_s}{\sqrt{\rho_f\left(2\gamma_{sf}\right)h_0}} \end{equation}

are the Ohnesorge numbers based on the film and surroundings viscosities, respectively. For this configuration, we keep ${Oh}_{f}$ constant at $0.05$ (based on the experiments of Reyssat & Quéré (Reference Reyssat and Quéré2006)), and vary the control parameter ${Oh}_{s}$ in § 5.

Figure 5. Computational domain for (a) two-phase and (b) three-phase Taylor–Culick retractions. For the classical case, panel (a) is used by replacing the surroundings (s) with air (a). The size of the domain is much larger than the hole radius $(\mathcal {L}_{max} \gg R(t))$. Furthermore, $\mathcal {L}_{max}/h_0 \gg \max ({Oh}_{f}, {Oh}_{s})$.

Note that the computational domain in figure 5(a) along with (4.5)–(4.8a,b) can be used to simulate classical Taylor–Culick retractions as well by replacing the surroundings (s) with air (a). We discuss the details of the classical configuration in Appendix A.

4.2.2. Three-phase Taylor–Culick configuration

In this configuration, we model the bursting of a water film at an oil drop-air interface by simulating the retraction of a fluid film ( f) on an initially flat oil bath (s), while ignoring the effects of the oil drop's curvature (as the retraction length in the early-time regime of figure 4b is much smaller than the oil drop radius, see figure 5b). We extend the traditional VoF method described in § 4.2.1 to tackle three fluids by using two VoF tracers: $\varPsi _1$, which is tagged as $1$ for the liquids (water film, f, and oil surroundings, $s$) and $0$ for air (a); and $\varPsi _2$ which is $1$ for the water film ( f) and $0$ everywhere else (figure 5b). Note that this implementation requires an implicit declaration of the surrounding phase (s), given by $\varPsi _2(1-\varPsi _1)$ (Sanjay et al. Reference Sanjay, Jain, Jalaal, van der Meer and Lohse2019; Sanjay Reference Sanjay2021a,Reference Sanjayb; Mou et al. Reference Mou, Zheng, Jian, Antonini, Josserand and Thoraval2021). Additionally, both $\varPsi _1$ and $\varPsi _2$ follow the VoF tracer advection equation,

(4.9)\begin{equation} \left(\frac{\partial}{\partial\tilde{t}} + \tilde{\boldsymbol{v}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}\right)\{\varPsi_1, \varPsi_2\} = 0, \end{equation}

and the dimensionless density ratio is (with $\rho _f = \rho _s$)

(4.10)\begin{equation} \tilde{\rho} = \varPsi_1 + \left(1-\varPsi_1\right)\left(\rho_a/\rho_f\right). \end{equation}

The Ohnesorge number ($Oh$) is now given by

(4.11)\begin{equation} Oh = \varPsi_1\varPsi_2{Oh}_{f} + \left(1-\varPsi_2\right)\varPsi_1{Oh}_{s} + \left(1-\varPsi_1\right){Oh}_{a}, \end{equation}

where ${Oh}_{f}$ and ${Oh}_{s}$ follow (4.8a,b), and ${Oh}_{a} = \eta _a/\sqrt {\rho _f(2\gamma _{sf})h_0}$ is the Ohnesorge number based on the viscosity of air. Both ${Oh}_{f}$ and ${Oh}_{a}$ are fixed at $10^{-1}$ and $10^{-3}$, respectively, for all the three-phase simulation data presented in this paper (see § 7), and we vary the control parameter ${Oh}_{s}$ in § 5. Lastly, the surface tension body force takes the form

(4.12)\begin{equation} \tilde{\boldsymbol{f}}_\gamma \approx \left(\gamma_{sa}/\gamma_{sf}\right)\left(\tilde{\kappa}_1/2\right)\tilde{\boldsymbol{\nabla}}\varPsi_1 + \left(\tilde{\kappa}_2/2\right)\tilde{\boldsymbol{\nabla}}\varPsi_2, \end{equation}

with $\gamma _{sa}$ and $\gamma _{sf}$ being the surface tension coefficients for the surroundings–air and surroundings–film interfaces, respectively.

Physically, such a configuration (figure 5b and (4.9)–(4.11)) ideally implies the presence of a zero thickness precursor film of the surrounding liquid (s, represented by $(1-\varPsi _2)\varPsi _1 = 1$, (4.11)) over the liquid film ( f, $\varPsi _1\varPsi _2 = 1$, (4.11)). Note that this numerical assumption is applicable only when it is thermodynamically favourable for one of the fluids (here s) to spread over all the other fluids, i.e. it has a positive spreading coefficient (de Gennes et al. Reference de Gennes, Brochard-Wyart and Quéré2004; Berthier & Brakke Reference Berthier and Brakke2012), $S \equiv \gamma _{af} - \gamma _{sf} - \gamma _{sa} > 0$, and the Neumann triangle collapses at the three-phase contact line. In reality, this precursor film will have a finite thickness controlled by microscopic forces (such as van der Waals forces; Vaynblat et al. (Reference Vaynblat, Lister and Witelski2001)), and is much smaller than the length scales that we can resolve numerically in the continuum framework. Indeed, for our numerical simulations, this precursor film has an effective thickness of $\varDelta /2$, where $\varDelta$ is the size of the finest grid employed in this work. We further assume that, on the time scale of film retraction, the effective spreading coefficient of the surrounding liquid (s) is $0$ (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). Consequently, the effective surface tension coefficient between the film and air is $\gamma _{af} = \gamma _{sf}+\gamma _{sa}$. This precursor film (Thoraval & Thoroddsen Reference Thoraval and Thoroddsen2013) is analogous to the mathematical model for spreading of a perfectly wetting liquid on a solid substrate (de Gennes et al. Reference de Gennes, Brochard-Wyart and Quéré2004; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), which regularizes the contact line singularity due to the numerical slip (with an effective slip length of $\Delta$/2; Afkhami et al. (Reference Afkhami, Buongiorno, Guion, Popinet, Saade, Scardovelli and Zaleski2018)) due to the discretization of the interface.

4.3. Note on non-dimensionalization in the viscous regime

For highly viscous surroundings (${Oh}_{s} > 1$), it is convenient to scale the velocities with the viscocapillary velocity scale $v_\eta$, owing to the dominant interplay between viscous and capillary stresses (Stone & Leal Reference Stone and Leal1989). Further, we can use the viscocapillary time $\tau _\eta$, film thickness $h_0$, and capillary pressure $p_\gamma$ to normalize the time, length and pressure dimensions, respectively,

(4.13a–c)\begin{equation} \tau_\eta = \frac{h_0}{v_\eta} = \frac{\eta_s h_0}{2\gamma_{sf}}, \quad v_\eta = \frac{2\gamma_{sf}}{\eta_s}, \quad p_\gamma = \frac{2\gamma_{sf}}{h_0}, \end{equation}

where $\gamma _{sf}$ is the surface tension coefficient between the film ( f) and the surroundings (s), $h_0$ the film thickness, and $\eta _s$ the viscosity of the surrounding medium. These viscocapillary scales modify the momentum equation as

(4.14)\begin{equation} \frac{\tilde{\rho}}{Oh_s^2}\left(\frac{\partial\tilde{\boldsymbol{v}}}{\partial \tilde{t}} + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left(\tilde{\boldsymbol{v}}\tilde{\boldsymbol{v}}\right)\right) ={-}\tilde{\boldsymbol{\nabla}} p + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left(2\tilde{\eta}\tilde{\boldsymbol{\mathcal{D}}}\right) + \tilde{\boldsymbol{f}}_\gamma. \end{equation}

Here, ${Oh}_{s}$ is the surroundings Ohnesorge number (4.8a,b), $\tilde {\rho }$ follows $\tilde {\rho } = 1$ and (4.10) for the two-phase and the three-phase configurations, respectively, and $\tilde {\boldsymbol {f}}_\gamma$ equals the corresponding expressions for the two configurations. Additionally, the dimensionless viscosities are given by

(4.15)\begin{equation} \tilde{\eta} = \begin{cases} \varPsi \left(\eta_f/\eta_s\right) + \left(1-\varPsi\right), & \text{two-phase case},\\ \varPsi_1\varPsi_2\left(\eta_f/\eta_s\right) + \left(1-\varPsi_2\right)\varPsi_1 + \left(1-\varPsi_1\right)\left(\eta_a/\eta_s\right), & \text{three-phase case}. \end{cases} \end{equation}

5. Taylor–Culick retractions: numerics

Figures 6 and 7 elucidate the two-phase and three-phase Taylor–Culick retractions. For low viscous surroundings (${Oh}_{s} \le 1$), figures 6(a) and 7(a) show the growth of the dimensionless hole radius ($\tilde {R}(t) = R(t)/h_0$) in time (normalized with the inertiocapillary time scale, $\tau _\gamma$), and the insets contain the growth rate of this hole: $\dot {\tilde {R}}_\gamma (t) = \tau _\gamma \,{\rm d}\tilde {R}(t)/{\rm d}t$. After the initial transients, the hole grows (i.e. the film retracts) linearly in time with a constant velocity ($v_f$). We can use this retraction velocity to calculate the film Weber number,

(5.1)\begin{equation} {We}_{f} \equiv \frac{\rho_fv_f^2h_0}{2\gamma_{sf}} = \lim_{\tilde{R} \to \infty}\dot{\tilde{R}}_\gamma^2, \end{equation}

which is represented by the black dashed lines in figures 6(a) and 7(a). Here ${We}_{f}$ is an output parameter of the retraction process. Note that, for very low ${Oh}_{s}$, as the rim grows with time, the inertial drag on the moving rim due to the surrounding medium overcomes the driving capillary forces, resulting in a decrease of the tip velocity (see insets of figure 6a and Jian et al. (Reference Jian, Deng and Thoraval2020b)). However, we can still calculate a velocity scale (and hence ${We}_{f}$) associated with the Taylor–Culick-like retraction immediately after the initial transients (as marked by the black dashed lines in the insets of figure 6a for the lowest ${Oh}_{s}$).

Figure 6. Two-phase Taylor–Culick retractions: temporal evolution of the dimensionless hole radius ($\tilde {R}(t)$) for (a) ${Oh}_{s} \le 1$ and (b) ${Oh}_{s} \ge 1$. Time is normalized using the inertiocapillary time scale, $\tau _\gamma = \sqrt {\rho _f h_0^3/\gamma _{sf}}$ in panel (a) and the viscocapillary time scale, $\tau _\eta = \eta _s h_0/\gamma _{sf}$ in panel (b). Insets of these panels show the variation of the dimensionless growth rate of the hole radius at different ${Oh}_{s}$, and mark the definitions of ${We}_{f}$ and ${Ca}_{s}$. Lastly, panel (c) illustrates the morphology of the flow at different ${Oh}_{s}$ at $\tilde {R} = 30$. In each snapshot, the left-hand side contour shows the velocity magnitude normalized with the (terminal) film velocity $v_f$ and the right-hand side shows the dimensionless rate of viscous dissipation per unit volume normalized using the inertiocapillary scales, represented on a $\log _{{10}}$ scale to differentiate the regions of maximum dissipation. Here, the film Ohnesorge number is ${Oh}_{f} = 0.05$. Also see supplementary movie 2.

Figure 7. Three-phase Taylor–Culick retractions: temporal evolution of the dimensionless hole radius ($\tilde {R}(t)$) for (a) ${Oh}_{s} \le 1$ and (b) ${Oh}_{s} \ge 1$. Time is normalized using the inertiocapillary time scale, $\tau _\gamma = \sqrt {\rho _f h_0^3/\gamma _{sf}}$ in panel (a) and the viscocapillary time scale, $\tau _\eta = \eta _s h_0/\gamma _{sf}$ in panel (b). Insets of these panels show the variation of the dimensionless growth rate of the hole radius at different ${Oh}_{s}$, and mark the definitions of ${We}_{f}$ and ${Ca}_{s}$. Lastly, panel (c) illustrates the morphology of the flow at different ${Oh}_{s}$ at $\tilde {R} = 30$. In each snapshot, the left-hand side contour shows the velocity magnitude normalized with the (terminal) film velocity $v_f$ and the right-hand side shows the dimensionless rate of viscous dissipation per unit volume normalized using the inertiocapillary scales, represented on a $\log _{{10}}$ scale to differentiate the regions of maximum dissipation. Here, the film Ohnesorge number is ${Oh}_{f} = 0.10$ and that of air is ${Oh}_{a} = 10^{-3}$. Also see supplementary movie 3.

Furthermore, when the surroundings is highly viscous (${Oh}_{s} \ge 1$), we plot the growing hole radius $\tilde {R}(t)$ as a function of time, which is normalized by the viscocapillary time scale $\tau _\eta$, see figures 6(b), 7(b) and § 4.3. The insets of these panels contain the growth rate of the hole, calculated as $\dot {\tilde {R}}_\eta = \tau _\eta \,{\rm d}\tilde {R}/{\rm d}t$. Once again, we observe that the growth of the hole (and the film retraction) depends linearly on time with a constant velocity, which can be used to calculate the surroundings capillary number,

(5.2)\begin{equation} {Ca}_{s} \equiv \frac{\eta_sv_s}{2\gamma_{sf}} = \lim_{\tilde{R} \to \infty}\dot{\tilde{R}}_\eta, \end{equation}

marked with the black dashed lines in figures 6(b), 7(b) and the corresponding insets. Here ${Ca}_{s}$ is another output parameter of the retraction process. Note that the velocity of the retracting film ($v_f$) is the same as the velocity scale in the surrounding medium ($v_s$), following the kinematic boundary condition at the circumference of the growing hole. Consequently, the two output parameters, ${We}_{f}$ (5.1) and ${Ca}_{s}$ (5.2), are related as ${Ca}_{s} = {Oh}_{s}\sqrt {{We}_{f}}$ (see § 7).

Lastly, figures 6(c) and 7(c) illustrate the flow morphologies for the two-phase and three-phase configurations, respectively, when the hole has grown to $\tilde {R} = 30$. Readers can refer to supplementary movies 2 and 3 for the temporal dynamics of the two-phase and three-phase configurations, respectively. Similar to the classical Taylor–Culick retraction case (Appendix A and supplementary movie 4), both the film and the surroundings move. However, unlike the classical case, even for low ${Oh}_{s}$, the surrounding medium takes away momentum from the film owing to inertia (added mass-like effect), thus reducing the retraction velocity (see insets of figures 6a and 7a). Furthermore, contrary to the classical case where the dissipation is highest at the neck connecting the rim to the rest of the film (see Appendix A), the dissipation in the other two configurations is spread out, and also occurs in the surrounding medium.

As the hole grows, the retracting film collects fluid parcels from upstream of the moving front and forms a rim (Culick Reference Culick1960; de Gennes Reference de Gennes1996; Villermaux Reference Villermaux2020). Essentially, the moving fluid parcels of the retracting tip collide with the fluid parcels upstream of the tip, which were initially at rest. The collisions are inelastic as a fraction of the available energy is dissipated by internal viscous fluid friction. For the classical and two-phase configurations, this rim entails a top–bottom symmetry, which is lost in the three-phase configuration. This is due to the air medium having significantly less inertia (added mass-like effect from the properties of the film) than the oil bath, causing the film to dig into the bath, and forming a hook-shaped rim (see figures 6ci–iii and 7ci–iii). Furthermore, the surrounding bath (s) engulfs the retracting film in order to feed the precursor film. This is a result of the high capillary pressure (high curvature) in the wedge region near the apparent three-phase contact line (Sanjay et al. Reference Sanjay, Jain, Jalaal, van der Meer and Lohse2019; Cuttle et al. Reference Cuttle, Thompson, Pihler-Puzović and Juel2021), which also aids in the formation of the hook-shaped rim (Peschka et al. Reference Peschka, Bommer, Jachalski, Seemann and Wagner2018). Moreover, for the cases where there is a density contrast between the film and the surroundings, the top–bottom symmetry can break down even for the two-phase configuration due to a flapping instability at very low ${Oh}_{s}$, as discussed by Lhuissier & Villermaux (Reference Lhuissier and Villermaux2009) and Jian et al. (Reference Jian, Deng and Thoraval2020b). Furthermore, as ${Oh}_{s}$ increases, the bulbous rim disappears, leading to slender, elongated retracting films. In the two-phase case, the retraction film maintains (top–bottom) symmetry (see figure 6ciii–v), and dissipation is highest in the viscous boundary layer in the surrounding medium (see figure 6c; further elaborated upon in § 8.2). However, the three-phase case features (top–bottom) asymmetric films owing to the accumulation of fluid towards the low-resistance air medium (see figure 7ciii–v), and dissipation is highest near the apparent three-phase contact line (see figure 7c; further elaborated upon in § 8.2). The disappearance of bulbous rim matches with the experimental observations (see § 3 and Reyssat & Quéré Reference Reyssat and Quéré2006).

Note that the numerical results presented in this section are complementary to the experiments on a film retracting in a submerged oil bath (two-phase case, Reyssat & Quéré (Reference Reyssat and Quéré2006)) and a film bursting at an air–liquid interface (three-phase case, § 3). The numerical simulations give us access to the cross-sectional view to elucidate the shape of the retracting films (figures 6 and 7), which is difficult to resolve experimentally. On the other hand, our axisymmetric (by definition) simulations do not show the azimuthal instabilities resulting in the corrugations at the oil–air–water contact line. Furthermore, as we focus only on the early-time dynamics, we also neglect the curvature of the oil drop in the case of the three-phase retractions. Nonetheless, we can still sufficiently compare the dependences of the retraction velocity on the Ohnesorge number ${Oh}_{s}$ of the surroundings (see § 7), along with the scaling relations that we develop in the next section for both the experimental and numerical datapoints.

6. Taylor–Culick retractions: a force perspective

The capillary and viscous forces, along with the inertia of the film and the surrounding media, govern the retraction dynamics. For the classical configuration (figure 1a), the viscosity and inertia of the outer medium are negligible. Furthermore, the film viscosity $\eta _f$ plays no role in determining the magnitude of the retraction velocity owing to the internal nature of the associated viscous stresses (Savva & Bush Reference Savva and Bush2009), as long as $Oh_{f}$ is less than the aspect ratio of the film (see Deka & Pierson Reference Deka and Pierson2020). Using these features, Taylor (Reference Taylor1959b) calculated $v_f = v_{{TC}}$ (1.1), resulting solely from momentum equilibrium, while disregarding the fate of the liquid accumulated in the rim (Villermaux Reference Villermaux2020). In terms of the dimensionless numbers introduced earlier (see § 5), (1.1) implies that ${We}_{f} = \rho _fv_f^2h_0/(2\gamma _{af})$ is constant and equal to 1 (see Appendix A for details of the retraction dynamics in the classical configuration). In this section, we delve into the different realizations of the dominating forces, and their implications, in the two-phase and three-phase configurations.

6.1. Two-phase Taylor–Culick retractions

For the two-phase configuration (figure 1b), if the viscosity of the oil phase is small (i.e. ${Oh}_{s} = \eta _s/\sqrt {\rho _f\gamma _{sf} h_0} \ll 1$), the Weber number based on the film velocity $v_f$ and the driving surface tension coefficient ($2\gamma _{sf}$), ${We}_{f} = \rho _fv_f^2h_0/(2\gamma _{sf})$ (5.1) has a value smaller than 1 (see inset of figure 6a). Nonetheless, the driving surface tension force $F_\gamma (t) \sim \gamma _{sf}(2{\rm \pi} R(t))$ (see figure 1b) still balances the inertial force $F_\rho (t) \sim \rho _f v_f^2(2{\rm \pi} R(t))h_0$. Note that since the oil (surrounding) and water (film) densities are very similar ($\rho _f \approx \rho _s$), we can still use $\rho _f$ for the density scale despite the added mass-like effect. Consequently, in this regime, the Weber number is still a constant during retraction (${We}_{f} \sim O(1)$, inset of figure 6a).

On the other hand, if the viscosity of the oil phase ($\eta _s$) is significantly higher (i.e. ${Oh}_{s} \gg 1$), the resistive viscous force $F_{\eta }(t)$ dominates over the inertial effects, as the surroundings Reynolds number $Re_s \equiv \rho _s v_s h_0/\eta _s \sim O(10^{-2}$). In such a scenario, the retraction dynamics will be governed by the balance between the capillary ($F_{\gamma }(t)$) and viscous ($F_{\eta }(t)$) forces (Fraaije & Cazabat Reference Fraaije and Cazabat1989; Reddy et al. Reference Reddy, Manivannan, Basavaraj and Thampi2020), given by

(6.1)\begin{equation} F_{\gamma}(t) \sim F_{\eta}(t), \end{equation}

where (from figure 1b)

(6.2)\begin{equation} F_\gamma(t) = 2 \gamma_{sf}\left(2{\rm \pi} R(t)\right). \end{equation}

For $F_{\eta }(t)$ in (6.1), one can consider the retracting rim to be a cylinder translating in a viscous flow (Reyssat & Quéré Reference Reyssat and Quéré2006; Eri & Okumura Reference Eri and Okumura2010). Thus, the viscous drag can be described by the Oseen approximation to the Stokes flow (Lamb Reference Lamb1975; Happel & Brenner Reference Happel and Brenner1983), which to the leading order is expressed as

(6.3)\begin{equation} F_{\eta}(t) \sim \eta_{s} v_f\left(2{\rm \pi} R(t)\right), \end{equation}

where the factor $2{\rm \pi} R(t)$ is due to the axisymmetric geometry. On equating (6.2) and (6.3), we get

(6.4)\begin{equation} v_f \sim \frac{\gamma_{sf}}{\eta_{s}} . \end{equation}

Moreover, $v_f = v_s$ (where $v_s$ is the velocity scale in the surrounding medium, see § 5). As a result, (6.4) implies that the capillary number ${Ca}_{s}$ (5.2) is constant, i.e.

(6.5)\begin{equation} {Ca}_{s} = \frac{\eta_s v_s}{2\gamma_{sf}} \sim O\left(1\right) . \end{equation}

Further, upon dividing both sides of (6.4) by the inertiocapillary velocity scale $v_\gamma = \sqrt {2\gamma _{sf}/(\rho _fh_0)}$ and squaring, we obtain

(6.6)\begin{equation} {We}_{f} \sim {Oh}_{s}^{{-}2} . \end{equation}

The aforementioned equations (6.5)–(6.6) denote the scaling laws for viscous two-phase Taylor–Culick retractions.

6.2. Three-phase Taylor–Culick retractions

For the three-phase configuration (figure 1c), in the viscous limit (${Oh}_{s} \gg 1$), the force balance is still given by (6.1), especially for the oils that are significantly more viscous than water. Here, the driving surface tension force can be expressed as (from figures 1c and 5b)

(6.7)\begin{equation} F_\gamma(t) = (\gamma_{sf} + \gamma_{sa} + \gamma_{sf} - \gamma_{sa}) \left(2{\rm \pi} R(t)\right)= 2 \gamma_{sf}\left(2{\rm \pi} R(t)\right), \end{equation}

assuming the presence of a precursor film of oil on top of the water film (see § 4.2.2, de Gennes et al. (Reference de Gennes, Brochard-Wyart and Quéré2004), Bonn et al. (Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009) and Thoraval & Thoroddsen (Reference Thoraval and Thoroddsen2013)). However, writing an expression for $F_{\eta }(t)$ is not as straightforward as the two-phase configuration. As can be observed from figure 7(c), during the retraction of the film, the oil climbs on top of the water, resulting in a strong flow in the wedge-like region close to the oil–air–water contact line. The rate of local viscous dissipation in this region is also very high (right-hand subpanels of figure 7c). Similar wedge flows have also been observed for moving contact lines on solid substrates (de Gennes Reference de Gennes1985; Marchand et al. Reference Marchand, Chan, Snoeijer and Andreotti2012; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). It has been reported that the wedge flow results in a viscosity-dependence of velocity that is weaker than 1/$\eta _{s}$ (Marchand et al. Reference Marchand, Chan, Snoeijer and Andreotti2012), but the exact nature of the dependence has hitherto not been quantified. The presence of a deformable liquid substrate on which the wedge flow occurs (the retracting water film in this case) complicates the situation even further – making it extremely difficult to arrive at the experimentally observed $v_f(\eta _{s}) \sim \eta _{s}^{-1/2}$ dependence (3.1) from a simple force balance. In § 8.2, we attempt to explain this scaling from an energetics point of view. Nevertheless, from the experiments, we know that the $v_f (\eta _{s})$ scaling is given by (3.1), which can be rewritten as

(6.8)\begin{equation} {Ca}_{s} \sim {Oh}_{s}^{1/2} . \end{equation}

Dividing both sides of (6.8) by $v_\gamma$ from (4.3a–c) and squaring, we obtain

(6.9)\begin{equation} {We}_{f} \sim {Oh}_{s}^{{-}1} . \end{equation}

Therefore, from (6.4), (6.6), (6.8) and (6.9), we hypothesize that the presence of the oil–air–water apparent contact line in the three-phase configuration dramatically alters the scaling relationships as compared with the two-phase configuration for ${Oh}_{s} > 1$ (see figures 6b and 7b). This will be further elaborated upon in § 8.2. Contrary to this scenario, for low ${Oh}_{s}$ numbers, the retraction velocities in both these configurations have the same scaling behaviour. Despite the presence of a hook-shaped rim in the three-phase case (figure 7ci-ii), we can still treat the moving rim and the surroundings as lumped elements. As a result, the driving surface tension force $\gamma _{sf}(2{\rm \pi} R(t))$ still balances the inertial force that scales with $\rho _f v_f^2(2{\rm \pi} R(t))h_0$, thus giving ${We}_{f} \sim O(1)$ (see figure 7a). In the next section, we demonstrate the validity of the scaling relations developed in this section.

7. Demonstration of the scaling relationships

Figure 8 illustrates the dependence of ${We}_{f}$ and ${Ca}_{s}$ on the Ohnesorge number ${Oh}_{s}$ of the surroundings for the retraction configurations described in figure 1. Note that the same datapoints are presented in both figures 8(a) and 8(b), following the relation ${Ca}_{s}~=~{Oh}_{s}\sqrt {{We}_{f}}$ (as $v_f = v_s$, see § 5). In figure 8(a), ${We}_{f} = 1$ marks the classical Taylor–Culick retraction limit, whereas for the two-phase and three-phase configurations, we identify two regimes: inertial (${Oh}_{s} < 1$) and viscous (${Oh}_{s} > 1$).

Figure 8. Regime maps visualized as ${We}_{f}$ versus ${Oh}_{s}$ in (a) and as ${Ca}_{s}$ versus ${Oh}_{s}$ in (b). The experimental datapoints (circles) correspond to the three-phase configuration (figure 1c) while the simulations (triangles) correspond to both the two-phase (figure 1b) and three-phase (figure 1c) configurations. The experimental datapoints (pentagrams) for the two-phase configuration have been adopted from Reyssat & Quéré (Reference Reyssat and Quéré2006) for their silicone oil–soap water (surroundings–film, s–f) dataset.

The inertial scaling is identical for both the two-phase and three-phase configurations: ${We}_{f} \sim O(1)$, which also implies ${Ca}_{s} \sim {Oh}_{s}$ (see §§ 6.1 and 6.2). The brown lines in figure 8 represent these two scaling relations.

The datapoints corresponding to the two-phase numerical simulations (from figure 6) are shown by the dark blue triangles. As ${Oh}_{s}$ increases, the retraction transitions from the inertial scaling (brown lines), to the viscous two-phase Taylor–Culick scaling: ${Ca}_{s}~\sim ~{Oh}_{s}^0$ (6.5) or ${We}_{f} \sim {Oh}_{s}^{-2}$ (6.6). We also plot the experimental datapoints from Reyssat & Quéré (Reference Reyssat and Quéré2006) for their silicone oil–soap water (surroundings–film, s–f) dataset, shown in figure 8 by the light blue pentagrams. In order to make these datapoints dimensionless, we use $h_0 = 100\,\mathrm {\mu }{\rm m}$ and $\gamma _{sf} = 7\,{\rm mN}\,{\rm m}^{-1}$, denoting the thickness of the soap film and the surroundings–film interfacial tension coefficient, respectively, in their experiments. We also neglect any Marangoni flow, or dynamic surface tension effects. Our simulations and scaling relationships are in reasonable agreement with the experimental datapoints of Reyssat & Quéré (Reference Reyssat and Quéré2006). Note that Reyssat & Quéré (Reference Reyssat and Quéré2006) tried to fit a trend line of $(\ln \eta _{s})/\eta _{s}$ (higher-order Oseen correction) through all of their experimental datapoints to obtain a good fit. When the same datapoints are plotted in figure 8(a,b), it is observed that some of their datapoints (corresponding to low ${Oh}_{s}$) are, in fact, in the transition between the inertial and the viscous regimes, while the rest of the datapoints show reasonable agreement with the viscous two-phase retraction dynamics given by (6.5) or (6.6).

We also plot the datapoints corresponding to our experiments (figure 4) and simulations (figure 7) for the three-phase configuration (figure 1c) in figure 8. We observe that, at low ${Oh}_{s}$, these datapoints follow the inertial dynamics (brown lines), while at higher ${Oh}_{s}$, the datapoints follow the scaling relationships given by (6.8) and (6.9): ${Ca}_{s} \sim {Oh}_{s}^{1/2}$ and ${We}_{f} \sim {Oh}_{s}^{-1}$, respectively (represented by the black lines in figure 8). Note that in order to non-dimensionalize the experimental datapoints shown in figure 4(c) (so that they can be plotted in figure 8), one needs to know the film thickness $h_0$. In the present experiments, the optical resolution was insufficient for accurate measurement of the film thickness prior to rupture. Moreover, as mentioned earlier, the breakup process itself is highly sensitive to experimental noise (see § 4.2 of Villermaux (Reference Villermaux2020)). Similar difficulties were also presumably experienced by Eri & Okumura (Reference Eri and Okumura2010) in their experiments of two-phase retraction, and they used a fitting parameter in their $v_f (\eta _{s})$ relation, which was a function of $h_0$. We also know from bubble bursting experiments (Doubliez Reference Doubliez1991; Lhuissier & Villermaux Reference Lhuissier and Villermaux2012) that the film thickness prior to breakup varies in the range $O(100\,{\rm nm})$–$O(10\,\mathrm {\mu }{\rm m})$. Moreover, in similar studies (Lhuissier & Villermaux Reference Lhuissier and Villermaux2012; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2012), the film thickness is retroactively calculated from the retraction velocity measurements. We can see from figure 8 that for low ${Oh}_{s}$, the dynamics are independent of the specific nature of the configuration (classical, two-phase or three-phase). Knowing $v_f$, $\eta _{s}$ and $\gamma _{sf}$, we can calculate the ${Ca}_{s}$ for the datapoint in figure 4(c) corresponding to $\eta _{s} = 4.94 \times 10^{-4}\,{\rm Pa}\,{\rm s}$. Fitting that ${Ca}_{s}$ value to the ${Ca}_{s} \sim {Oh}_{s}$ trend line (brown line) in figure 8(b), a value of $h_0 = 1.5\,\mathrm {\mu }{\rm m}$ can be calculated, which is within the range observed for previous experiments in a similar system (Lhuissier & Villermaux Reference Lhuissier and Villermaux2012). Using $h_0 = 1.5\,\mathrm {\mu }{\rm m}$ for the remaining experimental datapoints in figure 4(c) (for $\eta _{s} > 4 \times 10^{-3}$ Pa s) and calculating ${Ca}_{s}$, ${Oh}_{s}$ and ${We}_{f}$, we find that those datapoints (red circles) also collapse on the trend lines (black lines) along with the numerical simulations (yellow triangles) in figure 8. Note that a water film thickness of $h_0 = 1.5\,\mathrm {\mu }{\rm m}$ sets the ${Oh}_{f}$ at $0.1$ for the three-phase case, which is different from the ${Oh}_{f}$ that we use for the two-phase case (${Oh}_{f} = 0.05$ based on the experiments of Reyssat & Quéré (Reference Reyssat and Quéré2006)). Therefore, to justify comparison between the two cases, we varied ${Oh}_{f}$ in simulations from $0.01$ to $0.1$ and found that the dimensionless retraction velocities (${We}_{f}$ and ${Ca}_{s}$) are ${Oh}_{f}$-independent for both the two-phase and three-phase configurations (for ${Oh}_{f} < 1$). We also verify the ${Oh}_{f}$-independence experimentally by replacing the water in our bath by glycerol–water mixtures, and the measurements thus obtained (green circles) also follow the ${We}_{f} \sim {Oh}_{s}^{-1}$ and ${Ca}_{s} \sim {Oh}_{s}^{1/2}$ trendlines (black lines) in figures 8(a) and 8(b), respectively.

In summary, in §§ 6–7, we discussed the forces involved during the retraction of liquid films owing to the unbalanced capillary traction, followed by identification of the inertial (${Oh}_{s} < 1$) and viscous (${Oh}_{s} > 1$) regimes in the ${We}_{f}$ versus ${Oh}_{s}$ and ${Ca}_{s}$ versus ${Oh}_{s}$ dependences. We also checked the validity of the corresponding scaling behaviours in this section. To further understand the retraction dynamics due to the capillary traction, we focus on the different thermodynamically consistent energy transfer modes in the next section. Particularly, we try to understand the scaling relationship for the viscous three-phase Taylor–Culick retraction, that still eludes understanding from a momentum balance point of view (see § 8.2).

8. Taylor–Culick retractions: an energetics perspective

A retracting liquid sheet loses surface area and consequently releases energy (Dupré Reference Dupré1867, Reference Dupré1869; Rayleigh Reference Rayleigh1891; Culick Reference Culick1960), which further increases the kinetic energy of the system (i.e. film and surroundings). A part of this energy is lost in the process due to viscous dissipation. So, the overall energy budget entails

(8.1)\begin{equation} E_k^f(R) + E_k^s(R) + E_k^a(R) + E_\gamma(R) + E_d^f(R) + E_d^s(R) + E_d^a(R) = E_\gamma(R = 0), \end{equation}

where $E_\gamma$ is the surface energy, $E_k$ the kinetic energy and $E_d$ the viscous dissipation. The superscripts account for the film ( f), the surroundings (s) and air (a). Of course, for the two-phase case, the terms associated with air (a) do not exist as there is no air phase. Figure 9 depicts (8.1) for both the two-phase and three-phase configurations. Coincidentally, even for the three-phase configuration, the energies associated with the air phase are negligible (see figure 9, the dashed and dot–dashed lines overlap), even though the velocity field in air is not negligible (figure 7c). We keep $E^a_k(R)$ and $E^a_d(R)$ in the energy budget for the sake of completeness. In general, the sum of all these energies at any hole radius $R(t)$ equals the total surface energy at $R = 0$, i.e. the total energy available to the system. As the film retracts, it continuously releases energy, as its surface energy decreases. Therefore, to calculate (8.1), one can choose a reference for surface energy arbitrarily. In figure 9, the surface energy at a hole radius of $R = R_{{max}}$ is used as this arbitrary instance. This datum is chosen such that by the time the hole expands to $R_{{max}}$, the film would have reached a constant velocity. Furthermore, we can normalize the energies in (8.1) with the total surface energy released as the film retracts to a hole of radius $R_{{max}}$. The energy budget now reads

(8.2)\begin{equation} \bar{E}_k^f(\tilde{R}) + \bar{E}_k^s(\tilde{R}) + \bar{E}_k^a(\tilde{R}) + \Delta \bar{E}_\gamma(\tilde{R}) + \bar{E}_d^f(\tilde{R}) + \bar{E}_d^s(\tilde{R}) + \bar{E}_d^a(\tilde{R}) = 1. \end{equation}

Here, $\bar {E}(\tilde {R}) = E(\tilde {R})/(E_\gamma (0)-E_\gamma (\tilde {R}_{{max}}))$, $\Delta E_{\gamma }(\tilde {R}) = E_{\gamma }(\tilde {R}) - E_{\gamma }(\tilde {R}_{{max}})$, and $\tilde {R}=\tilde {R}(t)=R(t)/h_{0}$ is the dimensionless hole radius. The reader is referred to Appendix B for details of the energy budget calculations.

Figure 9. Energy budget at different ${Oh}_{s}$ for the (a) two-phase and (b) three-phase configurations. The energies $E$ are normalized by the total surface energy released as the film retracts creating a hole of radius $\tilde {R}_{{max}} = 100$ for $Oh_s \le 1$, and $\tilde {R}_{{max}} = 1000$ (two-phase case) and $\tilde {R}_{{max}} = 1200$ (three-phase case) for $Oh_s = 100$. Note that this $\tilde {R}_{{max}}$, and hence the surface energy datum, are arbitrarily chosen. We use hole radii that are large enough such that the sheets approach a constant velocity. The superscripts account for the film ( f), the surroundings (s) and air (a).

Removing the arbitrary datum described above and noting that there is a continuous injection of surface energy ($-\dot {E}_\gamma$, minus sign because the surface energy is decreasing with the growing hole) into the system, we can also write the energy budgets in terms of rates,

(8.3)\begin{equation} \dot{E}_k^f(\tilde{R}) + \dot{E}_k^s(\tilde{R}) + \dot{E}_d^f(\tilde{R}) + \dot{E}_d^s(\tilde{R}) = \left(-\dot{E}_\gamma(\tilde{R})\right). \end{equation}

Figure 10 visualizes (8.3) by plotting the proportion of the rate of surface energy released that goes into the rate of increase of kinetic energy and the rate of total viscous dissipation. From figures 10(ai–iii) and 10(bi–iii), we observe that these fractions saturate after initial transients. So, we also plot these steady state values (8.4) in figures 10(a) and 10(b) for the two-phase and three-phase configurations, respectively,

(8.4)\begin{equation} \left(\dot{E}/\dot{E}_{\gamma}\right)_\infty = \lim_{\tilde{R} \to \infty}\left(\frac{\dot{E}(\tilde{R})}{\dot{E}_{\gamma}(\tilde{R})}\right). \end{equation}

Figure 10. Variation of the rate of change of kinetic energy ($\dot {E}_k$) and viscous dissipation ($\dot {E}_d$) as proportions of the rate of energy injection ($-\dot {E}_\gamma$) with ${Oh}_{s}$ at steady state for the (a) two-phase and (b) three-phase configurations. For both cases, in the inertial limit (${Oh}_{s} \ll 1$), the fraction of energy that goes into kinetic energy and viscous dissipation are comparable. However, in the viscous limit (${Oh}_{s} \gg 1$), viscous dissipation in the surroundings dominates. Insets show the representative temporal variations of the ratio of the rate of change of energy ($\dot {E}$) to the rate of energy injection ($-\dot {E}_{\gamma }$) with dimensionless hole radius $\tilde {R}$ at three different ${Oh}_{s}$. The superscripts account for the film ( f), the surroundings (s) and air (a).

We devote the rest of this paper to understanding the distribution of the energy injection rate into the rates of increase of kinetic energy and viscous dissipation for both the inertial and viscous regimes.

8.1. Energy transfers in the inertial regime

We first focus on the energy balance in the classical Taylor–Culick retraction and the famous Dupré–Rayleigh paradox (Villermaux Reference Villermaux2020). Dupré (Reference Dupré1867, Reference Dupré1869) hypothesized that the total surface energy released during retraction manifests as the kinetic energy of the film (Rayleigh Reference Rayleigh1891). As a result, the predicted retraction velocity was off by a factor of $\sqrt {2}$ (see Appendix A), leading to discrepancies with experiments (Ranz Reference Ranz1959; Culick Reference Culick1960). Nonetheless, it is noteworthy that Dupré (Reference Dupré1867, Reference Dupré1869) reached the correct scaling relationship by identifying the essential governing parameters of classical sheet retractions.

Culick (Reference Culick1960) identified that the rate of surface energy released (8.9) should be distributed into an increase in kinetic energy of the rim and the viscous dissipation inside the film: $-\dot {E}_\gamma (t) = \dot {E}_k^f(t) + \dot {E}_d^f(t)$. The viscous dissipation can be attributed to the inelastic acceleration of the undisturbed film up to the velocity of the edge of the rim. Note that the dissipation is independent of the fluid viscosity and is given by (Culick Reference Culick1960)

(8.5)\begin{equation} \dot{E}_d^f(t) = \frac{1}{2}\frac{{\rm d}m(t)}{{\rm d}t}v_f^2 , \end{equation}

where $m(t)$ is the mass of the retracting film.

Coincidentally, this rate of viscous dissipation in the film is the same as the rate of increase in its kinetic energy (constant rim velocity; Culick (Reference Culick1960) and Villermaux (Reference Villermaux2020)). We confirm this hypothesis in Appendix A (see figures 12c,d, and Sünderhauf et al. (Reference Sünderhauf, Raszillier and Durst2002)), whereby

(8.6)\begin{equation} \dot{E}_k^f(t) \approx \dot{E}_d^f(t) \approx{-} \dot{E}_\gamma(t)/2. \end{equation}

Next, we delve into the energy transfers in the two-phase and three-phase configurations. In the inertial limit, in a manner akin to the classical case, the fraction of the rate of energy injection that goes into increasing the kinetic energy is similar to that of viscous dissipation. However, unlike the classical case, the kinetic energy as well as viscous dissipation are distributed among the film and the surrounding medium (figures 9 and 10, ${Oh}_{s} \ll 1$). We observe that

(8.7)\begin{equation} \left(\dot{E}_d^f(t) + \dot{E}_d^s(t)\right) \approx \left(\dot{E}_k^f(t) + \dot{E}_k^s(t)\right) \approx{-}\dot{E}_\gamma(t)/2. \end{equation}

In a manner reminiscent of Dupré (Reference Dupré1867, Reference Dupré1869), we can write

(8.8)\begin{equation} -\dot{E}_\gamma(t) \approx \left(\dot{E}_k^f(t) + \dot{E}_k^s(t)\right) \sim \left(\rho_f v_fh_0\left(2{\rm \pi} R(t)\right)\right)v_f^2, \end{equation}

where $v_f = v_s$ (kinematic boundary condition at the tip of the film) and $\rho _s = \rho _f$. Additionally, following Bohr & Scheichl (Reference Bohr and Scheichl2021) and Appendix B, the rate of change of surface energy is given by

(8.9)\begin{equation} \dot{E}_\gamma(t) \approx{-}F_{\gamma}(t) \frac{{\rm d}R(t)}{{\rm d}t} ={-}2\gamma_{sf}\left(2{\rm \pi} R(t)\right)v_f. \end{equation}

Using (8.8)–(8.9), and rearranging, we get

(8.10)\begin{equation} {We}_{f} = \frac{\rho_fv_f^2h_0}{2\gamma_{sf}} \sim O\left(1\right), \end{equation}

which is the same as the inertial scaling derived using the force balance (insets of figures 6a and 7a).

8.2. Demystifying dissipation in the viscous regime

In the viscous limit $({Oh}_{s} \gg 1)$, for both the two-phase and three-phase configurations, the surface energy released is entirely dissipated in the surrounding medium (figures 9 and 10), i.e.

(8.11)\begin{equation} -\dot{E}_\gamma(t) \sim \dot{E}_d^{s}(t). \end{equation}

In fact, this interplay between the surface energy and the viscous dissipation sets the velocity scale $(v_s)$ in the surrounding medium, which is equal to the retraction velocity ($v_f$, kinematic boundary condition at the hole). Therefore, to estimate this velocity, we first calculate the rate of viscous dissipation $\dot {E}_d^{s}(t)$, which depends on both the viscosity $\eta _{s}$ of the surrounding medium and the velocity gradients $\boldsymbol {\mathcal {D}}$, following the relation (see Appendix B)

(8.12)\begin{equation} \dot{E}_d^s = \int_{\varOmega_s} 2\eta_{s} \left(\boldsymbol{\mathcal{D}}:\boldsymbol{\mathcal{D}}\right) \,{\rm d}\varOmega_s = \int_{\varOmega_s}\varepsilon_s \,{\rm d}\varOmega_s. \end{equation}

Here, $\varepsilon _s$ is the rate of viscous dissipation per unit volume, and the integrals are evaluated over the volume $\varOmega _{s}$ of the surrounding medium. Note that $\varepsilon _s$ is the highest at the expanding hole, i.e. the tip of the retracting film in the case of two-phase retractions (figure 6c), and the macroscopic contact line in the case of three-phase retractions (figure 7c). The latter is analogous to wetting and dewetting of rigid surfaces (de Gennes Reference de Gennes1985; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). Motivated by this analogy, we calculate the local rate of viscous dissipation integrated over volume elements $\varOmega _s(r)$ centred at the expanding hole,

(8.13)\begin{equation} \dot{E}_d^{s, {local}}(r, t) = \int_{0}^{\varOmega_s(r)} \varepsilon_s(r, t) \,{\rm d}\varOmega_s, \end{equation}

where $r$ is the radial distance away from the hole (see insets of figure 11c). Additionally, in the viscous regime, we can use the viscocapillary velocity $v_\eta = 2\gamma _{sf}/\eta _s$ and the film thickness $h_0$ to non-dimensionalize (8.13) (see § 4.3 and Stone & Leal (Reference Stone and Leal1989)),

(8.14)\begin{equation} \dot{\tilde{E}}_d^{s, {local}}\left(\tilde{r}, \tilde{t}\right) \equiv \frac{\dot{E}_d^{s, {local}}\left(\tilde{r}, \tilde{t}\right)}{\eta_sv_\eta^2h_0} = \int_{0}^{\tilde{\varOmega}_s(\tilde{r})} \tilde{\varepsilon}_s\left(\tilde{r}, \tilde{t}\right) \,{\rm d}\tilde{\varOmega}_s. \end{equation}

Figure 11. Dissipation in the viscous limit (${Oh}_{s} \gg 1$) of Taylor–Culick retractions: evolution of the local rate of viscous dissipation $(\dot {\tilde {E}}_d^{s, {local}}(\tilde {r}, \tilde {t}))$ with dimensionless distance $\tilde {r} = r/h_0$ away from (a) the tip of the film in the two-phase configuration and (b) the macroscopic three-phase contact line in the three-phase configuration. In insets (ii), this distance is normalized with the dimensionless viscous boundary layer thickness in the surrounding medium, $\tilde {\delta }_\nu = \delta _\nu /h_0 = {Oh}_{s}\sqrt {\tilde {t}}$. Here, $\tilde {R} = R/h_0$ and $\tilde {t} = t/\tau _\eta$ are the dimensionless hole radius and dimensionless time, respectively. (c) Variation of the total viscous dissipation rate per unit circumference of the hole $(\dot {\tilde {E}}_d^{s}(\tilde {t})/(2{\rm \pi} \tilde {R}(\tilde {t})))$ at steady state with the surroundings Ohnesorge number ${Oh}_{s}$.

Figures 11(ai) and 11(bi) show that the local viscous dissipation increases as we move away from the hole (increasing $\tilde {r}$). Furthermore, the energy dissipated increases in time as the region of flow expands, owing to the increasing hole radius and the dominant radial flow. To rationalize this increase, we plot the rate of local viscous dissipation per unit circumference of the hole in figures 11(aii) and 11(bii).

For the two-phase case, the viscous dissipation occurs in the viscous boundary layer $( \tilde {\delta }_\nu \sim {Oh}_{s}\sqrt {\tilde {t}})$ and saturates at $\tilde {r} \approx \tilde {\delta }_\nu$ (figure 11aii). However, for the three-phase case, we can identify two distinct regions of viscous dissipation, the wedge region close to the macroscopic contact line, where the viscous dissipation per unit circumference of the expanding hole increases steeply $(\tilde {r} < 0.01\tilde {\delta }_\nu )$, and the viscous boundary layer $(\tilde {r} < 0.1\tilde {\delta }_\nu )$, beyond which it saturates (figure 11bii). Furthermore, this saturation value gives the total viscous dissipation per unit circumference of the hole,

(8.15)\begin{equation} \frac{\dot{\tilde{E}}_d^s(\tilde{t})}{\left(2 {\rm \pi}\tilde{R}(\tilde{t})\right)} = \lim_{\tilde{r} \to \infty}\frac{\dot{\tilde{E}}_d^{s, {local}}\left(\tilde{r}, \tilde{t}\right)}{\left(2 {\rm \pi}\tilde{R}(\tilde{t})\right)}, \end{equation}

which is shown in figure 11(c) as a function of ${Oh}_{s}$. We observe that for the two-phase case, the total dissipation is independent of ${Oh}_{s}$, whereas in the three-phase case, it scales with ${Oh}_{s}^{1/2}$:

(8.16)\begin{equation} \dot{\tilde{E}}_d^{s}(\tilde{t}) \sim \begin{cases} {Oh}_{s}^{0}\left(2{\rm \pi}\tilde{R}(\tilde{t})\right), & \text{two-phase case},\\ {Oh}_{s}^{1/2}\left(2{\rm \pi}\tilde{R}(\tilde{t})\right), & \text{three-phase case}. \end{cases} \end{equation}

Moreover, upon non-dimensionalizing (8.9) using the same scales as used in (8.14), and noting that $v_f = v_s$ and ${Ca}_{s} = \eta _sv_s/(2\gamma _{sf})$, we get

(8.17)\begin{equation} -\dot{\tilde{E}}_\gamma(t) \equiv \frac{\gamma_{sf} v_f\left(2 {\rm \pi}R(t)\right)}{\eta_sv_\eta^2h_0} = {Ca}_{s} \left(2 {\rm \pi}\tilde{R}(\tilde{t})\right). \end{equation}

Lastly, equating (8.16) and (8.17), we get

(8.18)\begin{equation} Ca_s \sim \begin{cases} {Oh}_{s}^{0}, & \text{two-phase case},\\ {Oh}_{s}^{1/2}, & \text{three-phase case}. \end{cases} \end{equation}

In summary, in this section, we confirmed our hypothesis that the presence of the oil–air–water contact line in the three-phase configuration dramatically alters the scaling relationships and dynamics as compared with the two-phase configuration (see § 6.2). We also relate the dimensionless retraction velocity ${Ca}_{s}$ with the control parameter ${Oh}_{s}$ in the viscous limit by following the location and magnitude of the local rate of viscous dissipation during Taylor–Culick retractions in viscous surroundings.