Taylor–Culick retractions and the influence of the surroundings

Abstract

When a freely suspended liquid film ruptures, it retracts spontaneously under the action of surface tension. If the film is surrounded by air, the retraction velocity is known to approach the constant Taylor–Culick velocity. However, when surrounded by an external viscous medium, the dissipation within that medium dictates the magnitude of the retraction velocity. In the present work, we study the retraction of a liquid (water) film in a viscous oil ambient (two-phase Taylor–Culick retractions), and that sandwiched between air and a viscous oil (three-phase Taylor–Culick retractions). In the latter case, the experimentally measured retraction velocity is observed to have a weaker dependence on the viscosity of the oil phase as compared with the configuration where the water film is surrounded completely by oil. Numerical simulations indicate that this weaker dependence arises from the localization of viscous dissipation near the three-phase contact line. The speed of retraction only depends on the viscosity of the surrounding medium and not on that of the film. From the experiments and the numerical simulations, we reveal unprecedented regimes for the scaling of the Weber number ${We}_{f}$ of the film (based on its retraction velocity) or the capillary number ${Ca}_{s}$ of the surroundings versus the Ohnesorge number ${Oh}_{s}$ of the surroundings in the regime of large viscosity of the surroundings (${Oh}_{s} \gg 1$), namely ${We}_{f} \sim {Oh}_{s}^{-2}$ and ${Ca}_{s} \sim {Oh}_{s}^{0}$ for the two-phase Taylor–Culick configuration, and ${We}_{f} \sim {Oh}_{s}^{-1}$ and ${Ca}_{s} \sim {Oh}_{s}^{1/2}$ for the three-phase Taylor–Culick configuration.

1. Introduction

Liquid films, sheets and shells have piqued the interest of fluid dynamicists for nearly two centuries (Savart 1833a,b,c; Taylor 1959a,b; Clanet 2007; Villermaux 2020). A freely suspended liquid film warrants further attention among these configurations as it is inherently metastable owing to its high surface area. Indeed, if a large enough (Taylor & Michael 1973; Villermaux 2020) hole nucleates on the film, the sheet will spontaneously retract to reduce its surface area. Such interfacial destabilization leading to film rupture and bursting can also result in waterborne disease transmission (Bourouiba 2021). The bursting of liquid films at an oil–air interface is important for various industrial applications in the chemical and petrochemical engineering sectors as well. One area of particular interest is underwater oil spills in oceans, such as the Deepwater Horizon spill in 2010 in the Gulf of Mexico (Summerhayes 2011). For these spills, droplets (or slugs) of oil may rise to the free surface of water via buoyancy, and then rupture the free surface of water directly above it. The water film will retract upon rupture, and the oil will spread on the water surface, thus perpetuating an environmental hazard. Perhaps the most widely studied example of sheet destabilization and retraction is during the bursting of liquid (e.g. soap) films in air – an area of active research since the pioneering works (Dupré 1867, 1869; Rayleigh 1891; Taylor 1959b; Culick 1960; McEntee & Mysels 1969) in the late nineteenth and mid-twentieth century, to the more recent investigations (Bremond & Villermaux 2005; Müller, Kornek & Stannarius 2007; Lhuissier & Villermaux 2012; Munro et al. 2015; Deka & Pierson 2020). In these studies, the outer medium is assumed passive (inviscid and zero-inertia). The origin of the nucleation of the initial hole in the film can be manifold (Lohse & Villermaux 2020). After film rupture, the internal viscous stresses in the film do not contribute to the momentum balance, but dictate the distribution of momentum within the film (Savva & Bush 2009), as long as the Ohnesorge number of the film (ratio of its viscocapillary to inertiocapillary time scales, see § 4) is less than its aspect ratio (see Deka & Pierson 2020). Nonetheless, half of the surface energy released goes into internal viscous dissipation (see Appendix A and Culick 1960; de Gennes 1996; Sünderhauf, Raszillier & Durst 2002; Villermaux 2020).

A representative schematic of the situation mentioned above is shown in figure 1(a) (henceforth referred to as the classical Taylor–Culick configuration), where the water film ( f) of thickness $h_0$ is retracting in air (a) under the action of surface tension. The retraction velocity, $v_f$, in such a scenario is constant (after a period of initial transience) and approaches the Taylor–Culick velocity given by

(1.1)\begin{equation} v_{{TC}} = \sqrt{\frac{2\gamma_{af}}{\rho_f h_0}} , \end{equation}

where $2\gamma _{af}$ is the net surface tension driving the retraction ($\gamma _{af}$ being the interfacial tension coefficient between the film and air) and $\rho _f$ is the density of the liquid film (figure 1a). Furthermore, it was observed that during the retraction, the liquid collects in a thicker rim at the retracting edge of the film, particularly for low viscosity liquids (Rayleigh (1891), Ranz (1959), Pandit & Davidson (1990); not depicted in figure 1b). The seminal work of Keller (1983) further explores the retraction of these films of non-uniform thickness.

Figure 1. Schematics depicting the configurations studied in the present work: (a) retraction of a water film ( f) of thickness $h_0$ in an air (a) environment (classical configuration); (b) retraction of water film ( f) in an oil (s) environment (two-phase configuration); (c) retraction of a water film sandwiched between air and oil (three-phase configuration). The dot–dashed line represents the axis of rotational symmetry, and $R(t)$ is the radius of the growing hole centred at this axis. In all the schematics, the water film is retracting from left to right with velocity $v_f$, as indicated by the arrow, and $\gamma _{ij}$ denotes the surface tension coefficient between fluids $i$ and $j$.

The effect of viscosity of the film ($\eta _f$) during its retraction process has also been studied (Debrégeas, Martin & Brochard-Wyart 1995; Debrégeas, de Gennes & Brochard-Wyart 1998). Brenner & Gueyffier (1999) showed that although viscosity does not have any effect on the constant retraction velocity, it can have a significant effect on the shape of the retracting edge of a planar film. They report that if the radial extent of the film is greater than its Stokes length ($= \eta _f/(\rho _f v_{{TC}})$), a growing rim is formed, whereas the rim is absent for the converse situation. Savva & Bush (2009) extended the work by Brenner & Gueyffier (1999) for highly viscous films, and also developed a lubrication model for the retraction dynamics of a circular hole. They concluded that although viscosity does not determine the magnitude of the constant retraction velocity, it does dictate the time required (post rupture) to attain that constant velocity, which increases with increasing viscosity. Recently, Pierson et al. (2020) and Deka & Pierson (2020) revisited the viscous retraction dynamics by exploring self-similar solutions for slender filaments and sheets of finite length.

The rheological properties of the film also influence the retraction dynamics (Dalnoki-Veress et al. 1999; Tammaro et al. 2018; Villone, Hulsen & Maffettone 2019; Kamat, Anthony & Basaran 2020). For instance, Sen et al. (2021) showed that viscoelastic filaments can retract at velocities higher than the Newtonian Taylor–Culick limit due to elastic tension. Moreover, the retraction dynamics of liquids have also been studied in the context of dewetting for a wide range of scenarios (Redon, Brochard-Wyart & Rondelez 1991; Brochard-Wyart, Martin & Redon 1993; Shull & Karis 1994; Andrieu, Sykes & Brochard 1996; Lambooy et al. 1996; Haidara, Vonna & Schultz 1998; Buguin, Vovelle & Brochard-Wyart 1999; Péron, Brochard-Wyart & Duval 2012; Peschka et al. 2018; Kim et al. 2020). Lastly, there has also been a recent surge in the study of liquid retraction in other configurations, such as liquid strips (Lv, Clanet & Quéré 2015), smectic films (Trittel et al. 2013), foam films (Petit, Le Merrer & Biance 2015) and emulsion films (Vernay, Ramos & Ligoure 2015).

In all the aforementioned studies, the surrounding medium is assumed to play no role in the rupture dynamics. A question naturally arises: What happens when the outer medium also interacts with the retracting film? In particular, how do the viscosity and the inertia of the outer medium influence the rupture dynamics (Mysels & Vijayendran 1973; Joanny & de Gennes 1987; Reyssat & Quéré 2006; Jian, Deng & Thoraval 2020b)? A representative schematic for such a scenario is shown in figure 1(b), where a water film is retracting in a viscous oil ambient. This geometry will henceforth be referred to as the two-phase configuration where the net surface tension force responsible for retraction is $2\gamma _{sf}$ ($\gamma _{sf}$ being the interfacial tension coefficient between the film and the surrounding medium, see figure 1b). In such a situation, viscous dissipation is not limited only to the retracting film, but is also present in the ambient. If the ambient happens to be significantly more viscous than the film, then the dissipation in the ambient dominates. In such a situation, the retraction velocity is still a constant. However, unlike the classical case, the velocity depends on the viscosity $\eta _{s}$ of the ambient medium (Martin, Buguin & Brochard-Wyart 1994; Reyssat & Quéré 2006). Common realizations of this configuration include relaxation of filaments and droplets in a viscous medium (Stone & Leal 1989), or that of air-films during drop impact (Jian et al. 2020a,b). Additionally, in this context, Anthony, Harris & Basaran (2020) showed that the so-called ‘inertially limited viscous regime’ in the early-times of drop coalescence (Paulsen et al. 2012; Paulsen 2013) stems from a Taylor–Culick-type retraction of the air film between the deformable drops.

In the present work, we study the influence of the surrounding medium on the retraction velocity of a ruptured liquid film using both force balance and energy conservation arguments. To accomplish this goal, along with the two canonical configurations shown in figures 1(a) and 1(b), we also study the retraction dynamics of a liquid film sandwiched between air and a viscous oil bath. A representative schematic is shown in figure 1(c). This geometry will henceforth be referred to as the three-phase configuration. This paper elucidates this case experimentally by inflating an oil droplet at the water–air interface and letting the water film rupture. Such a configuration can also be found in the early stages of water film retraction when an air bubble approaches a water–oil interface if the oil layer is thick enough (Feng et al. 2014, 2016). Furthermore, we also use direct numerical simulations to demystify the retraction dynamics by using a precursor film-based three-fluid volume of fluid (VoF) method. We show that the film in this three-phase configuration still retracts with a constant velocity, and similar to the two-phase case, the retraction velocity depends on the viscosity $\eta _{s}$ of the oil bath. However, this dependence is weaker in the three-phase configuration. Furthermore, we reveal an unprecedented scaling relationship for the retraction velocity of the film, which arises from the localization of the viscous dissipation near the three-phase contact line.

The paper is organized as follows. Section 2 describes the problem statement for the three-phase Taylor–Culick retractions along with the experimental method employed to probe this configuration. The results from these experiments are discussed in § 3. Section 4 presents the numerical framework, and § 5 describes the simulation results for both the two-phase and three-phase configurations. Section 6 demonstrates the balance of forces in Taylor–Culick retractions, followed by the corresponding scaling relationships in § 7. Further, § 8 analyses the overall energy balance, highlighting the differences in the viscous dissipation mechanisms between the two-phase and three-phase configurations. The work culminates with conclusions in § 9. Throughout the manuscript, we refer to Appendix A for discussions on the classical Taylor–Culick retractions, and use the experimental datapoints from Reyssat & Quéré (2006) for the two-phase configuration.

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 2013).

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 2001), 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. (2016). 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 2012) 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 (2020)). 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 2012; de Maleprade, Clanet & Quéré 2016). The rupture location may also be determined by a ‘prehole’ formation, also observed by Vernay et al. (2015) for the bursting of emulsion-based liquid sheets. In their work, Vernay et al. (2015) 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 (2020). 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 1990; Brenner & Gueyffier 1999; Sünderhauf et al. 2002). 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. 2020), 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é 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é 2006; Oldenziel et al. 2012) 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 2008; Kavehpour 2015). The nature of the corrugations is reminiscent of the sharp tips observed during selective withdrawal (Cohen & Nagel 2002; Courrech du Pont & Eggers 2006, 2020) or tip streaming (Montanero & Gañán Calvo 2020). 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é (2006) for the instabilities observed in their experiments for film retraction in the two-phase configuration. Tseng & Prosperetti (2015) 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 (1999) 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. 1994; Eri & Okumura 2010). 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 2013–2022; Sanjay 2021b), which uses the one-fluid approximation (Tryggvason, Scardovelli & Zaleski 2011) 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 1992).

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. 1992)

(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 2009). We follow the same sign convention as Tryggvason et al. (2011, 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é (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. 2019; Sanjay 2021a,b; Mou et al. 2021). 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. 2004; Berthier & Brakke 2012), $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. (2001)), 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. 2009). Consequently, the effective surface tension coefficient between the film and air is $\gamma _{af} = \gamma _{sf}+\gamma _{sa}$. This precursor film (Thoraval & Thoroddsen 2013) is analogous to the mathematical model for spreading of a perfectly wetting liquid on a solid substrate (de Gennes et al. 2004; Bonn et al. 2009), which regularizes the contact line singularity due to the numerical slip (with an effective slip length of $\Delta$/2; Afkhami et al. (2018)) 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 1989). 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}

4.4. Domain size and boundary conditions

Figure 5 depicts the computational domains. The left-hand boundary represents the axis of symmetry with origin marked at $(0, 0)$. We set no-penetration and free-slip boundary conditions to all other domain boundaries along with zero gradient conditions for pressure. These boundaries are far away from the expanding hole and do not affect its growth. Furthermore, the size of the domain is chosen such that $\mathcal {L}_{{max}} \gg \text {max}({Oh}_{f}, {Oh}_{s})$, with a minimum $\mathcal {L}_{{max}}$ of $200$ for ${Oh}_{s} \ll 1$. We have varied this domain size to ensure that the simulations are independent of its value. Note that, if this condition is not met, the assumption of infinite film, which is essential for the theoretical scaling relations developed in this work, will fail (Deka & Pierson 2020).

We employ adaptive mesh refinement to correctly resolve the different interfaces as well as regions of high velocity gradients (and hence, high viscous dissipation, see Appendix B). To ensure that the velocity field is captured accurately, these refinement criteria (see Sanjay 2021b) effectively maintain a minimum of $40$ cells across the thickness of the film (i.e. $h_0/\Delta \ge 40$). As the apparent three-phase contact line and the viscous boundary layer are critical in the present work, the refinement criteria maintain a minimum of $40$ cells in the wedge region near the apparent three-phase contact line. Furthermore, the viscous boundary layer is almost $10$ times larger than the film thickness (see § 8.2). Consequently, a minimum of $400$ cells in the viscous boundary layer in the surrounding medium is needed to properly resolve the velocity gradients. We have conducted extensive grid independence studies so that the final results (energy transfers and the retraction velocity) are independent of the number of grid cells.

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. (2020b)). 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 1960; de Gennes 1996; Villermaux 2020). 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. 2019; Cuttle et al. 2021), which also aids in the formation of the hook-shaped rim (Peschka et al. 2018). 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 (2009) and Jian et al. (2020b). 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é 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é (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 2009), as long as $Oh_{f}$ is less than the aspect ratio of the film (see Deka & Pierson 2020). Using these features, Taylor (1959b) calculated $v_f = v_{{TC}}$ (1.1), resulting solely from momentum equilibrium, while disregarding the fate of the liquid accumulated in the rim (Villermaux 2020). 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 1989; Reddy et al. 2020), 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é 2006; Eri & Okumura 2010). Thus, the viscous drag can be described by the Oseen approximation to the Stokes flow (Lamb 1975; Happel & Brenner 1983), 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. (2004), Bonn et al. (2009) and Thoraval & Thoroddsen (2013)). 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 1985; Marchand et al. 2012; Snoeijer & Andreotti 2013). 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. 2012), 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é (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é (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é (2006). Note that Reyssat & 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 (2020)). Similar difficulties were also presumably experienced by Eri & Okumura (2010) 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 1991; Lhuissier & Villermaux 2012) 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 2012; Thoroddsen et al. 2012), 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 2012). 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é (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é 1867, 1869; Rayleigh 1891; Culick 1960), 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 2020). Dupré (1867, 1869) hypothesized that the total surface energy released during retraction manifests as the kinetic energy of the film (Rayleigh 1891). As a result, the predicted retraction velocity was off by a factor of $\sqrt {2}$ (see Appendix A), leading to discrepancies with experiments (Ranz 1959; Culick 1960). Nonetheless, it is noteworthy that Dupré (1867, 1869) reached the correct scaling relationship by identifying the essential governing parameters of classical sheet retractions.

Culick (1960) 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 1960)

(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 (1960) and Villermaux (2020)). We confirm this hypothesis in Appendix A (see figures 12c,d, and Sünderhauf et al. (2002)), 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é (1867, 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 (2021) 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 1985; Bonn et al. 2009; Snoeijer & Andreotti 2013). 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 (1989)),

(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.

9. Conclusion and outlook

In this paper, we have studied the effects of the surrounding media on the retraction dynamics of liquid sheets in three canonical configurations. In the classical Taylor–Culick configuration, the interplay between capillarity and inertia of the film results in a constant retraction velocity. We can further neglect the surrounding medium as it does not influence the retraction process. However, for a film retracting in a dense and viscous oil (two-phase configuration), and that at an oil–air interface (three-phase), both inertia and viscosity of the oil phase influence the retraction process. The former presents itself as an added mass-like effect. Even though capillarity still governs the constant retraction velocity, the surrounding medium's inertia reduces the magnitude of the film's momentum as it retracts.

Moreover, when the viscosity of the oil is significantly higher than that of the film, the viscous stresses dictate the retraction process and set the velocity scale. To further demystify the energy balance in this process, we have used thermodynamically consistent energy transfer mechanisms to understand the fate of the released surface energy owing to the loss of surface area of the retracting film. This energy is injected into the system and manifests itself as kinetic energy and viscous dissipation. In the inertial regime, the proportions of kinetic energy and viscous dissipation are the same, conforming to the analyses of Culick (1960). However, in the viscous regime, the total surface energy released goes into viscous dissipation in the surroundings.

Following the lumped elements analysis, motivated by Taylor (1959b) and Culick (1960), we also developed scaling relations to relate the non-dimensionalized retraction velocity (${We}_{f}$ and ${Ca}_{s}$) with the control parameter ${Oh}_{s}$. In the inertial limit, the Weber number ${We}_{f}$ based on the retraction velocity is a constant for all three configurations. On the other hand, in the viscous limit, the retraction velocity in the two-phase configuration scales with the viscocapillary velocity scale (constant capillary number, ${Ca}_{s} \sim O(1)$); while for the three-phase configuration, the capillary number ${Ca}_{s}$ increases with increasing ${Oh}_{s}$, owing to the localization of viscous dissipation near the three-phase contact line.

A natural extension of the present work would be to understand the retraction of non-Newtonian sheets and filaments (Sen et al. 2021) in similar surroundings. In such scenarios, the retraction dynamics will depend not only on capillarity and viscosity as described in this work, but also on the rheological properties of both the film and the surroundings. Furthermore, in a broader perspective, the precursor film-based three-fluid VoF method can be used to elucidate several spreading phenomena, both at small and large scales, e.g. drop-film interactions in the inkjet printing process (Lohse 2022) and late time spreading during oil spillage (Hoult 1972), respectively.