2.1. Theoretical framework and numerical modelling

The one-fluid formulation of incompressible Navier–Stokes equations (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999) is employed to model the two-phase flow represented by the liquid curtain interacting with the initially quiescent gaseous environment:

(2.1a)$$\begin{gather} \dfrac{\partial u^\star_i}{\partial x^\star_i}=0, \end{gather}$$

(2.1b)$$\begin{gather}\rho\left( \dfrac{\partial u^\star_i}{\partial t^\star} + u^\star_j\,\dfrac{\partial u^\star_i}{\partial x^\star_j}\right)=-\dfrac{\partial p^\star}{\partial x^\star_i} + \rho g_i + \dfrac{\partial}{\partial x^\star_j} \left[\mu\left(\dfrac{\partial u^\star_i}{\partial x^\star_j} +\dfrac{\partial u^\star_j}{\partial x^\star_i} \right)\right] + F_\sigma, \end{gather}$$

(2.1c)$$\begin{gather}\dfrac{\partial C}{\partial t^\star} + \dfrac{\partial C u^\star_i}{\partial x^\star_i}=0. \end{gather}$$

The vectors $\boldsymbol {u}^\star = (u^\star,v^\star,w^\star )$ and $\boldsymbol {g} = (g,0,0)$ represent the flow velocity and the gravity acceleration, respectively, while $p^\star$ is the pressure. The force per unit volume due to the surface tension is denoted as $F_\sigma = \sigma \kappa ^\star n_i \delta _S$, where $\sigma$ is the surface tension coefficient, $\kappa ^\star$ is the mean gas–liquid interface curvature, and $\boldsymbol {n} = (n_x,n_y,n_z)$ is the outward pointing normal vector to the interface. The Dirac distribution function $\delta _S$ is equal to 1 at the interface, and 0 elsewhere. Note that all the dimensional quantities, except the gravity acceleration $g$ and the fluid material properties ($\rho _l$, $\rho _a$, $\mu _l$, $\mu _a$, $\sigma$), are denoted with the superscript $\star$. The density $\rho$ and viscosity $\mu$ fields are discontinuous across the interface separating the two fluids, namely

(2.2a)$$\begin{gather} \rho=\rho_a + (\rho_l -\rho_a)C , \end{gather}$$

(2.2b)$$\begin{gather}\mu=\mu_a + (\mu_l -\mu_a)C , \end{gather}$$

where subscripts $l$ and $a$ refer to liquid and ambient phases, respectively, and the volume fraction field $C$ is equal to either $1$ or $0$ in the liquid and gaseous regions.

The system (2.1a)–(2.1c) is solved using the finite volume method in the open-source code BASILISK (http://basilisk.fr), an improved version of Gerris (Popinet Reference Popinet2003) that has been used extensively and validated for several liquid jet flow problems (among others, see Della Pia, Chiatto & de Luca Reference Della Pia, Chiatto and de Luca2020; Schmidt & Oberleithner Reference Schmidt and Oberleithner2020; Agbaglah Reference Agbaglah2021b; Schmidt et al. Reference Schmidt, Tammisola, Lesshafft and Oberleithner2021). The code employs the volume-of-fluid method by Hirt & Nichols (Reference Hirt and Nichols1981) to track the interface on an octree structured grid, allowing for adaptive mesh refinement based on a criterion of wavelet-estimated discretization error (van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018), with no special treatment required in the presence of liquid phase breakup. The calculation of the surface tension term is based on the balanced continuum surface force technique (Francois et al. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006), which is coupled with a height-function curvature estimation to avoid the generation of spurious currents. For exhaustive details about BASILISK, the reader is referred to Popinet (Reference Popinet2003, Reference Popinet2009).

2.2. Geometrical layout and problem formulation

A schematic description of the computational domain is reported in figure 1(a), where the gravity direction is denoted by $x$. The domain is a cubic region extending along the $x$ (streamwise), $y$ (transverse) and $z$ (spanwise) directions, where $L^\star = 50 H^\star _i$ is the domain side length. The spatial coordinates $x$, $y$ and $z$ have been scaled with respect to the inlet (i.e. at the streamwise station $x=0$) sheet thickness $H^\star _i$ (namely $x=x^\star /H^\star _i$), while the corresponding streamwise $u$, transverse $v$ and spanwise $w$ velocity components have been made dimensionless with respect to the inlet mean liquid velocity $U^\star _i$ ($u = u^\star /U^\star _i$).

Figure 1. (a) Schematic representation of the computational domain, and (b) inflow velocity profile (2.3a). The gravity $g$ is directed along the streamwise direction $x$; $y$ is the lateral coordinate, and $z$ is the spanwise one.

The curtain issues into an initially quiescent gaseous environment (blue region in figure 1) from a rectangular slot of dimensions $H^\star _i \times W^\star _i$, representing the initial thickness (along the transverse $y$ coordinate) and width (along the spanwise $z$ direction) of the sheet, respectively. The origin of the reference frame coincides with the centre of the slot, and the curtain shape is initialized as a parallelepiped with volume $L^\star \times H^\star _i \times W^\star _i$ (red region in figure 1a), which is employed to start the computation. Dirichlet boundary conditions are enforced at the inlet: following Kacem (Reference Kacem2017), in the liquid region ($-1/2< y<1/2$, $-{{AR}}/2< z<{{AR}}/2$, where ${{AR}} = W^\star _i/H_i^\star$ is the sheet aspect ratio), a fully developed parabolic velocity profile with a proper error function modification is imposed to account for viscous effects at the slot boundaries ($y = \pm 1/2$, $z = \pm {{AR}}/2$), and the conditions read

(2.3a)$$\begin{gather} u=\dfrac{3}{2}\left(1 - 4y^2 \right) \text{erf} \left(\dfrac{{{AR}}}{2} - z\right) \text{erf} \left(\dfrac{{{AR}}}{2} + z\right), \end{gather}$$

(2.3b)$$\begin{gather}v = w=0, \end{gather}$$

(2.3c)$$\begin{gather}C=1 . \end{gather}$$

On the remaining part of the inlet plane, namely within the gaseous phase, a no-slip condition is imposed; the inflow velocity profile $u(x=0,y,z)$ is represented in figure 1(b). The four lateral boundary planes ($y= \pm 25$, $z= \pm 25$) are equipped with homogeneous Neumann boundary conditions for all variables, and on the outlet plane ($x=50$) a standard outflow condition

(2.4a)$$\begin{gather} \dfrac{\partial u}{\partial x} = \dfrac{\partial v}{\partial x}=\dfrac{\partial w}{\partial x} = \dfrac{\partial C}{\partial x}=0, \end{gather}$$

(2.4b)$$\begin{gather}p=0, \end{gather}$$

is considered. The same $u(y,z)$ profile enforced as inlet boundary condition (2.3a) is also employed as the initial velocity distribution throughout the entire sheet length. We note here explicitly that a slight modification of the initial and boundary conditions presented above is required to obtain the last results reported in § 4.3. In particular, to simulate an edge-free curtain flow, homogeneous Neumann boundary conditions are replaced by periodic conditions on the lateral walls of the domain. Moreover, the parallelepiped curtain shape employed as initial condition is extended all along the spanwise direction, thus matching the whole domain side length.

The computational domain is discretized with an adaptive mesh up to a maximum number of $2^9$ grid points along each spatial dimension, corresponding to a minimum mesh size $\varDelta ^\star \approx H^\star _i/11$, and approximately $134$ million cells if a uniform grid was used; the mesh resolution effect on the results hereafter presented is reported in Appendix A. Note that recently, a resolution $\varDelta ^\star / H^\star _i \approx 6$ has been shown by Agbaglah (Reference Agbaglah2021a) to be valuable in capturing holes expansion and collision in a thin liquid sheet.

The $n=10$ physical quantities involved in the problem ($\rho _l$, $\rho _a$, $\mu _l$, $\mu _a$, $g$, $U^\star _i$, $W^\star _i$, $H^\star_i$, $L^\star$, $\sigma$) can be rearranged in $m=n-3$ independent dimensionless numbers: a possible choice is reported in table 1 ($r_\rho$, $r_\mu$, $\varepsilon$, ${{AR}}$, $Re$, $Fr$, $We$), which also specifies a set of numerical values corresponding to the reference configuration analysed in § 3. Following Chubb et al. (Reference Chubb, Calfo, McConley, McMaster and Afjeh1994), the Froude number is based on the curtain width $W_i^\star$. Note that a different arrangement of the dimensional parameters leads to the definition of the Ohnesorge number $Oh = \mu /\sqrt {2 H^\star _i \rho _l \sigma }$, whose value is also reported in table 1.

Table 1. Dimensionless governing parameters and their values for the reference case.

4.1. Supercritical regime

The unsteady behaviour of the curtain shape in the $xz$ plane influenced by the time evolution of the hole is shown qualitatively in figure 6 and quantified in figure 7. The temporal evolution of the curtain shape is depicted in figure 6 in terms of volume fraction $C$ contour maps. Following the approach by Karim et al. (Reference Karim, Suszynski and Pujari2021), the trajectories of three representative points of the hole are monitored in figure 7, namely the northernmost $x_N(t)$ (red), southernmost $x_S(t)$ (blue) and centre $x_{C}(t)$ (green) points of the hole. The analysis is performed for a relatively low supercritical Weber number value $We=1.1$, for which surface tension is expected to have a relatively strong influence on the hole-driven curtain dynamics. The aspect ratio is ${{AR}}=40$. Note that the temporal variable $t$ in figures 6 and 7 is scaled with respect to the reference time $t^\star _{ref} = H^\star _i/u^\star _c$, the hole being introduced at $t=0$ in the steady curtain configuration.

Figure 6. Maps of volume fraction $C$ in the $xz$ plane at different time instants. The liquid phase is represented in red, the gaseous phase in blue. Here, $We = 1.1$, ${{AR}}=40$, with (a) $t = 0$, (b) $t = 1$, (c) $t = 2.2$, (d) $t = 3.3$, (e) $t = 4.9$, (f) $t = 6.7$, (g) $t = 7.4$, (h) $t = 8.6$, (i) $t = 10.2$.

Figure 7. (a) Hole initial shape in the $xz$ plane. (b) Curtain shape around the hole region in the $xy$ plane. (c) Hole northernmost $x_N(t)$ (red curve), southernmost $x_S(t)$ (blue curve) and centre $x_{C}(t)$ (green curve) points trajectories in the $xz$ plane. In (a,b), the liquid phase is shaded in red, with the gaseous phase in blue. Here, $We = 1.1$, ${{AR}} = 40$.

Following the sequence in figure 6, it can be noted that the initially circular hole is advected downstream progressively while undergoing a continuous shape deformation, until it closes and leaves the sheet intact. In the early stages of the evolution ($t < 3.3$), the temporal variation of the curtain shape is influenced only by the competition between the hole rim retraction, due to the surface tension acting on the hole perimeter, and the hole fall motion, determined mainly by the streamwise-oriented gravity force. As for the $x_N(t)$ motion, the surface tension counteracts gravity, the two forces having opposite directions in the northernmost point of the hole, while they are both directed downstream in the southernmost one (see figures 7a,b). This determines a progressive stretching of the hole along the streamwise direction (figures 6a–c), with a falling velocity $\dot{x}_N = 1.0$ smaller than $\dot{x}_S = 3.6$, as can be appreciated by the different initial slopes of the red and blue curves in figure 7(c); the streamwise velocity of the hole centre assumes a value $\dot{x}_{C} = 2.2$ (slope of the green curve in figure 7c). It is straightforward to verify that the orders of magnitude of these velocities are consistent with each other. As a matter of fact, if one assumes that the velocity of the hole centre point is essentially the falling velocity of the main liquid stream, then the velocities of the northernmost and southernmost points of the hole can be estimated by subtracting and adding, respectively, to the centre point velocity the ‘net’ rim retraction velocity, whose order of magnitude is given by the Taylor–Culick model (which is of unit order in dimensionless scale).

At intermediate stages ($3.3< t<8.6$), the hole dynamics starts to be influenced by the lateral edges of the curtain, namely the hole no longer expands freely. The curtain edges slow down the overall expansion (figures 6d–g), determining a reduction of the velocity $\dot{x}_S$ from 3.6 to 0.9 (figure 7(c) at $t=3.3$). The hole thus assumes a triangular shape, analogously to the underlying liquid sheet. Finally, during the last stages of the evolution ($t > 8.6$), the hole expansion is counteracted not only by the curtain edges, but also by the approaching vertical tail (figures 6h,i). Therefore, the streamwise stretching and the overall expansion of the hole stop: the southernmost point velocity becomes negative ($\dot{x}_S = -2.3$; see figure 7(c) at $t=9.0$), thus leading to the closure of the hole before it reaches the outlet section of the physical domain, i.e. the curtain length.

The hole dynamics in $yz$ planes is shown qualitatively in figure 8 and quantified in figure 9. Figures 8(a,b) depict the temporal evolution of the hole rim shape of the right hemi-curtain (i.e. for $z>0$) in two $yz$ planes located at the streamwise stations $x=13$ and $x=19$, respectively. Note that the geometrical features of the shapes significantly recall the elliptic, cigar-like and peanut-like shapes found in the pioneering combined theoretical–experimental investigation by Chubb et al. (Reference Chubb, McConley, McMaster and Afjen1993). The spanwise positions of the hole rim tip $z_T$ (leftmost point of the hole, denoted by the arrows in figure 8) and the retraction velocity $w_T= \dot{z}_T$ at the same streamwise stations are shown in figures 9(a,b) and 9(c,d), respectively, for $x=13$ and $x=19$, where the retraction velocity is scaled with respect to the Taylor–Culick reference value, i.e. $w_T = w^\star _T/u^\star _c$. Note that the time interval considered in each panel is computed as $\Delta t = t-t^0_h(x)$, where $t^0_h(x)$ is the time instant when the hole reaches the streamwise $x$ location. Since the hole grows progressively, moving along the streamwise direction for $x \in [0,20]$, as shown in figures 6 and 7 discussed previously, the time interval $\Delta t$ at $x=19$ (figures 8b and 9c,d) is greater than the corresponding interval at $x=13$ (figures 8a and 9a,b).

Figure 8. Temporal evolution of the hole rim shape in the $yz$ planes for the right hemi-curtain (i.e. $z>0$) for (a) $x=13$, and (b) $x=19$. The liquid phase is shaded in red, the gaseous phase in blue. The arrows denote the hole tip. Here, $We=1.1$, ${{AR}} = 40$.

Figure 9. Spanwise position $z_T$ and velocity $w_T$ of the hole rim tip in $yz$ planes (right hemi-curtain for $z>0$) as a function of time, for (a,b) $x=13$, and (c,d) $x=19$. The dashed lines in (b,d) denote the values $w^\star _T/u^\star _c= \pm 1$. Here, $We=1.1, {{AR}} = 40$.

The time evolution of the sheet shape reported in figures 8(a,b) exhibits a fast retraction of the rim at the early instants, at both $x=13$ and $x=19$, after its abrupt formation when the hole reaches the two streamwise stations. This can be appreciated qualitatively by looking at the large distance between the green and red curves at $y=0$, and it is quantified by the relatively higher slope of curves $z_T$ and $w_T$ at small $\Delta t$ values in figure 9. At later stages, the rim motion slows down, and the tip position $z_T$ reaches a maximum (for $\Delta t = 1.2$ and $2.1$, at $x=13$ and $19$, respectively) and then decreases. This corresponds to an inversion of the retraction process, due to the advection of the hole along the streamwise direction (blue and black curves in figure 8), which finally leaves both the $yz$ planes at relatively large time instants ($z_T = 0$ for $\Delta t = 2.5$ and $5.5$, at $x=13$ and $19$, respectively).

The main difference between the hole curtain dynamics at the two streamwise stations considered here is due to the increasing effect of the curtain edges as $x$ increases. At the upstream station, $x=13$, the edges are relatively far from the hole, which expands freely under the effect of surface tension. The hole-driven rim dynamics is thus characterized by a monotonic decreasing trend of the retraction velocity $w_T$, which eventually assumes negative values when the hole progressively leaves the station. On the contrary, the hole-induced rims are relatively close to the curtain edges at the downstream station, $x=19$, causing a mutual interaction that produces an oscillatory behaviour of the velocity $w_T$ in the range $w_T \in [-1,1]$, before the hole finally approaches the curtain tail.

To justify the relatively high values of the retraction velocity $w_T$ of the tip point, one should consider that the hole is immersed in a flow field subjected to the gravity acceleration. As a consequence, the ‘total’ retraction velocity of the left tip of the hole is given by the substantial derivative

(4.2)\begin{equation} w_T = \dfrac{\partial z_T}{\partial t} + u\,\dfrac{\partial z_T}{\partial x}, \end{equation}

where the unsteady term ${\partial z_T}/{\partial t}$ represents the pure retraction contribution usually modelled by the Taylor–Culick velocity, and the convective term $u ({\partial z_T}/{\partial x})$ is the tip displacement rate due to the local slope of the hole (in the $xz$ plane) advected downstream with velocity $u$. The convective term is indeed significant only when the hole tip at the given station $x$ corresponds to the maximum of the absolute value $|\partial z_T /\partial x|$ (i.e. the hole tip in the $yz$ plane coincides with the northernmost hole point in the $xz$ plane). On the contrary, the convective derivative is negligible when the hole tip coincides with the hole leftmost point in the $xz$ plane, so that the order of magnitude of $w_T$ coincides with the Taylor–Culick velocity, namely $w_T = \partial z_T/\partial t = {O}(1)$.

To sum up, the main result arising from the analysis performed in this section is that in supercritical conditions, the hole perturbation introduced in the curtain does not influence the dynamics of the flow at large times. Indeed, whilst the hole undergoes an initial transient expansion, its interaction with the curtain borders determines first a slowdown and then a reversal of the expansion process, thereby producing the hole closure upstream of the sheet exit section (see also supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.543).

4.2. Transcritical regime

The unsteady dynamics of the curtain in transcritical conditions is studied by considering two values of the Weber number close to the critical unit threshold, namely $We=1.1$ (as in § 4.1) and $We=0.95$.

The time evolution of the hole area $A^\star$ with respect to its initial value $A^\star _i$ is reported in figure 10, respectively for $We=1.1$ (black curves) and $We=0.95$ (red curves). The trend of both curves reported in figure 10(a) is the same: the hole area first increases due to the hole expansion at the early stages of its evolution, reaches a maximum at the peak time $t=t_{p}$, and then decreases to zero due to the hole closure process. Although the supercritical and subcritical behaviours are qualitatively analogous, two main differences can be detected. First, both the peak time $t_{p}$ and the corresponding area amplification $A^\star _p/A^\star _i$ are greater for $We=1.1$ ($t_{p} = 3.3$, $A^\star _p/A^\star _i = 58$) than for $We=0.95$ ($t_{p} = 2.7$, $A^\star _p/A^\star _i = 44$). Moreover, the initial amplification rate is higher for $We=0.95$ rather than $We=1.1$, as can be appreciated more clearly by looking at figure 10(b). The two curves are characterized by an exponential growth at the early time instants, $A^\star /A^\star _i = {\rm e}^{\alpha t}$, with the hole expansion rate $\alpha$ increasing from $\alpha =5/2$ to $\alpha =4$ by reducing the Weber number from supercritical to subcritical conditions. It is interesting to note that an exponential growth of hole amplifications inside two-dimensional liquid sheets was also found in the combined theoretical–numerical investigation by Savva & Bush (Reference Savva and Bush2009), who performed simulations in the high viscous limit ($Oh \gg 1$) and did not consider the effect of gravity on the hole dynamics. On the other hand, here $Oh \ll 1$ (see table 1 in § 2.2) and gravitational effects are fully taken into account.

Figure 10. (a) Time evolution of the hole area $A^\star$ with respect to its initial value $A^\star _i$ at the supercritical-to-subcritical flow transition; (b) zoom around the first time instants. Black curves denote $We=1.1$, and red curves denote $We=0.95$. The peak time $t_p$ and the corresponding area amplification are also highlighted in (a) for the case $We=1.1$. Here, ${{AR}}=40$.

Analogies and differences between the $We > 1$ and $We < 1$ cases can be explained easily by recalling the discussion in § 4.1 about the physical mechanisms involved in the hole dynamics. At the early stages of the evolution, the hole expands freely under the effect of the surface tension acting on its rims, which is clearly stronger in the subcritical rather than supercritical regime, thus determining a faster expansion of the area in the former case. On the other hand, as soon as the hole approaches the central vertical region of the curtain, the surface tension acting on the curtain edges counteracts the hole expansion, the effect being again stronger for $We < 1$ than for $We > 1$, thus determining an earlier decay (and consequently a lower maximum of the hole area) of the subcritical curve. The same considerations can be retrieved by looking at figure 11, which reports the hole northernmost ($x_N(t)$) and southernmost ($x_S(t)$) points trajectories for $We=1.1$ and $We=0.95$ (black and red curves, respectively). Moreover, the analysis of figure 11 reveals that in the transition from supercritical to subcritical flow regimes, the hole northernmost point velocity (slope of $x_N(t)$ curves) decreases slightly. Scaled with respect to the Taylor–Culick reference value $u^\star _c$, the velocity reads as $\dot{x}^\star _N/u^\star _c = 1.0$ for $We=1.1$, and $\dot{x}^\star _N/u^\star _c =0.8$ for $We=0.95$, respectively.

Figure 11. Hole northernmost $x_N(t)$ and southernmost $x_S(t)$ points trajectories in supercritical-to-subcritical flow transition conditions. Here, ${{AR}}=40$.

4.3. Edge-free subcritical regime

In the previous sections we observed that the dynamics of a hole introduced in the falling liquid sheet is influenced crucially by the presence of the lateral free edges of the curtain itself, determining a slowdown of the hole advected by the main stream, in both supercritical and weakly subcritical regimes. In this respect, there is no discontinuity in the sheet behaviour when traversing the critical Weber threshold. On the other hand, it has also been found that for finite values of the aspect ratio ${{AR}}$, the steady curtain flow configuration can no longer be achieved for lower $We$ numbers (e.g. for $We<0.8$), due to the spontaneous appearance of various irregular holes interacting with each other inside the sheet. This occurrence will be illustrated in Appendix B.

To address the analysis of the free dynamics of the hole, namely not influenced by the lateral edges of the sheet, the configuration of infinite ${{AR}}$ is investigated in the present section, for various subcritical $We$ numbers. This study has been performed by locating the hole perturbation (4.1) at a station further upstream, namely $x_h=2$. Moreover, by employing periodic boundary conditions on the lateral walls of the domain (see the discussion reported in § 2.2), the edge-free curtain flow has been simulated.

The time evolution of the curtain shape in the $xz$ plane is reported in figure 12 for the slightly subcritical Weber number value equal to $We=0.9$. It can be seen that in the absence of the curtain borders, the hole expands in a wide region (figures 12a–c) while it is convected downstream by the gravitational acceleration (figures 12d–f). As time increases further, the hole is expelled progressively at the outlet section (figures 12g,h), leaving the liquid curtain in its original unperturbed configuration (figure 12i). It is also worth noticing that by simulating the edge-free curtain flow in supercritical conditions ($We=1.1$), we retrieved the same overall behaviour outlined in figure 12. For the sake of brevity, edge-free supercritical curtain flow simulations have not been reported herein.

Figure 12. Time evolution of the volume fraction $C$ in the $xz$ plane of an edge-free curtain in slightly subcritical conditions ($We=0.9$). Time $t$ increases from (a) to (i), as (a) $t = 0$, (b) $t=5.8$, (c) $t = 9.8$, (d) $t = 11.8$, (e) $t = 14.8$, (f) $t = 17.8$, (g) $t = 21.9$, (h) $t = 31.9$, (i) $t = 39.9$. The liquid phase is represented in red, the gaseous one in blue.

For the same edge-free condition, the investigation has been extended to various lower values of subcritical Weber number. The results are summarized in figure 13 in terms of the hole northernmost point trajectory $x_N(t)$ as a function of $We$. Note that the initial streamwise position of the point is $x_N(0)=1.5$ for all the cases considered, being the hole perturbation located at $x_h=2$. By inspection of figure 13, it can be seen that after an initial transient period approximately equal to $t=1.5$, all curves assume a remarkably linear trend. This time $t=1.5$ corresponds to 98 time steps of the temporal numerical integration, and it represents approximately the time needed by the computation to fully develop the toroidal interface curvature around the hole. Note also that due to the low Ohnesorge number of the present case, $Oh=0.002$, following Savva & Bush (Reference Savva and Bush2009) and Pierson et al. (Reference Pierson, Magnaudet, Soares and Popinet2020), one can admit that this time interval is sufficient for the rim retraction velocity to reach the asymptotic value of the Taylor–Culick model.

Figure 13. Hole northernmost point trajectory $x_N(t)$ by decreasing the Weber number in subcritical edge-free conditions ($We <1$). For each $We$, the black dashed line represents the least squares linear fitting of the corresponding numerical data for $t>1.5$.

The major result is that as $We$ decreases, the velocity $\dot{x}_N$ (slope of the dashed curves in figure 13) decreases progressively until it changes sign, becoming negative for $We < 0.4$ (see also figure 14, where the velocity is scaled with respect to the Taylor–Culick reference value $u^\star _c$). Note that we verified the domain-size independency of results reported in figure 14, focusing attention around the threshold value $We=0.4$. We performed simulations (not reported herein) for both smaller ($L=40$) and larger ($L=60$) domains than the one considered in this work, and no appreciable variations of the trend $x_N(We)$ in the range $We \in [0.3,0.5]$ were detected. The validity of the results of figure 14 can be verified by considering that the figure reports the velocity of the rim top point, which is due to the gravitational velocity minus the Taylor–Culick one (in dimensionless terms, this last being of unit order). As a consequence, for instance, for $We_1=0.2$ and $We_2=0.9$, the gravitational velocities are, respectively, equal to $u_1=0.5+1=1.5$ and $u_2=-0.3+1=0.7$. Since the Weber number scales with the square of the gravitational velocity, equality must hold between $We_2/We_1= 4.5$ and $u_2^2 / u_1^2= 4.59$. Below the threshold value $We=0.4$, the surface tension force becomes strong enough to reverse the falling motion of the hole, which retracts upstream towards the domain inlet while expanding along the $z$ (lateral) direction. This is shown clearly in supplementary movie 2, which reports the entire hole expansion process for $We=0.2$. The upward retraction of the hole finally determines the breakup of the curtain, resulting in the columnar shape shown in figure 15, which resembles the various classic experimental findings of the literature (among others, see Pritchard Reference Pritchard1986; de Luca & Meola Reference de Luca and Meola1995; Marston et al. Reference Marston, Thoroddsen, Thompson, Blyth, Henry and Uddin2014) and the more recent contribution by Kyotoh et al. (Reference Kyotoh, Sekine and Roknujjaman2022). We also verified that the liquid columns pattern shown in figure 15 is unaffected qualitatively by increasing the domain size up to $L=60$ (not reported). For example, the spanwise wavelength of the pattern, which is approximatively $\lambda _z \approx 10$, does not change by increasing the domain size.

Figure 14. Hole northernmost point velocity (scaled with respect to the Taylor–Culick reference value $u^\star _c$) as a function of the Weber number in subcritical edge-free conditions ($We < 1$). The red dashed line denotes the zero value.

Figure 15. Time evolution of the edge-free three-dimensional curtain shape (frontal view in the $xz$ plane) in breakup conditions ($We=0.2$; see also supplementary movie 2). The liquid phase is represented in white, the gaseous phase in blue. Time $t$ increases from (a) to (d).

It is evident that the above discussion conceptually recovers the celebrated stability criterion of Brown (Reference Brown1961), namely that the condition of inlet $We > 1$ prevents the sheet rupture, although this criterion is very raw because it does not take into account the local competition between surface tension, driving the rim retraction, and the incoming momentum flux at the hole top point. Sunderhauf et al. (Reference Sunderhauf, Raszillier and Durst2002) extended the validity of Brown's criterion by studying in detail the hole retraction dynamics. In particular, they found that the hole top point first moves downwards and reaches a maximum vertical displacement at a certain instant; after this instant, the top margin of the hole starts to move upwards and initiates the rupture. However, gravity effects were not included in their analysis. Furthermore, Sunderhauf et al. (Reference Sunderhauf, Raszillier and Durst2002) stated that the rupture is avoided when, before the ‘turning’ instant, the hole top point is swept out of the curtain flow domain. This conclusion is precisely what we have found in this investigation. Moreover, the gravity effects have been taken into account fully here, also including the influence of the curtain aspect ratio up to the limit case of edge-free sheet flow. As a final remark, we note that it will be interesting to investigate the effect of gravity on the hole-induced rim transverse corrugations, similar to those leading to the formation of liquid fingers (Agbaglah, Josserand & Zaleski Reference Agbaglah, Josserand and Zaleski2013; Agbaglah & Deegan Reference Agbaglah and Deegan2014). Such transverse instabilities could indeed develop along the streamwise (gravity) direction with a large enough computational domain.