On the coupling instability of a gas jet impinging on a liquid film

Abstract

We investigate the dynamics of a gas jet impinging perpendicular to a thin liquid film dragged by a rising vertical substrate. This configuration is relevant to the jet-wiping process in hot-dip galvanization and it is unstable. Previous studies analysed the dynamics of the instability in the case of liquids with low Kapitza numbers (highly viscous liquids), more amenable to experimental and numerical investigations. This work extends the previous investigations by focusing on the wiping at much higher Kapitza numbers, which are more relevant to the galvanizing process. The simulations are carried out by combining volume of fluid and large-eddy simulations, and the dynamics of the gas–liquid interaction is analysed using extended multiscale proper orthogonal decomposition. The simulations allowed for analysing the jet-wiping instability in new flow conditions. Despite the largely different conditions, the results show that the interaction between the gas jet and the liquid film is qualitatively similar, featuring two-dimensional waves in the liquid correlated with oscillations and deflections of the gas jet in all cases. The wave characteristics (e.g. frequency and propagation speed) scale remarkably well using the Shkadov-like scaling based on the liquid, suggesting a dominant role of the liquid film in the coupling, and potentially enabling extrapolation of the results to a broader range of wiping conditions. Finally, we use the numerical results to discuss the limitations of liquid-film models, which constitute currently the only possible approach to study the jet-wiping process in industrial conditions.

1. Introduction

Thin liquid films interacting with gas flows are frequently encountered in engineering applications and in nature. The dynamics of these systems have a significant impact on the performance of industrial equipment such as evaporators, condensers and chemical reactors (Hanratty & Engen 1957; McCready & Hanratty 1985), as well as in processes like coating (Buchlin 1997) or arc welding (Berghmans 1972). Using direct numerical simulations, Lombardi, De Angelis & Banerjee (1996) and De Angelis, Lombardi & Banerjee (1997) have shown that, in some configurations, the large difference of scales in the dynamics of the two flows allows for treating the liquid phase as a ‘rigid’ boundary for the gas phase, thus decoupling their dynamics. However, this assumption is only valid when the interfacial deformations (waves) are confined to the viscous sublayer of the gas (Rosskamp, Willmann & Wittig 1998). When this condition is not met, as in most configurations of practical interest, understanding the gas–liquid interaction is crucial to improving the performance of the two-phase system. Examples of configurations with strongly coupled gas-interface interaction are the cavity flow produced by a jet impinging on a liquid bath (Ojiako et al. 2020) or the waves in falling liquid films subject to concurrent or counter-current laminar gas flows (Dietze & Ruyer-Quil 2013).

A third example, which is the focus of this work, is the jet-wiping process. It consists in using impinging planar gas jets to control the thickness of a thin film dragged by an upward-moving substrate. The process is widely used in hot-dip galvanizing because it allows for a precise and contactless control of the amount of zinc deposited on a steel substrate for corrosion protection. Simplified analytical models to predict the average coating thickness downstream jet wiping have been proposed by Thornton & Graff (1976), Buchlin (1997) and Yoneda (1993), while early stability analyses based on the Orr–Sommerfeld problem (Tuck 1983), and integral models (Tu & Ellen 1986) have concluded that the configuration is neutrally stable. This result was later confirmed by Hocking et al. (2011) using quasisteady models also accounting for the gas-induced shear stress at the interface, and by Ivanova et al. (2023) and Barreiro-Villaverde, Gosset & Mendez (2021) who investigated the dynamics of waves on thin films dragged by a solid wall. Nevertheless, all these works assumed a stationary jet flow or no jet flow. Experimental and numerical investigations by Gosset (2007) and Myrillas et al. (2013), later extended by Gosset, Mendez & Buchlin (2019) and Mendez, Gosset & Buchlin (2019b), have shown that the jet wiping is intrinsically unstable and characterized by a strong coupling between the gas and liquid phases, resulting in large oscillations of the impinging jet combined with large waves in the liquid film. This instability produces non-uniformities (long-wavelength undulations) on the final coating, limiting the quality of products requiring high surface finish, such as those for the automotive industry.

Besides the industrial relevance, understanding the origin of the jet-wiping instability remains a fundamental and unsolved problem in fluid dynamics, especially because this flow configuration remains inaccessible to comprehensive experimental and numerical investigations in industrial conditions. The link between the unsteadiness of the jet flow and the response of the coating layer was analysed via integral boundary layer (IBL) models introduced by Kapitza (1948) and Shkadov (1970) in a quasisteady formalism by Johnstone et al. (2019) and using self-similar and weighted integral boundary layer (WIBL) models (Ruyer-Quil & Manneville 2000, 2002) by Mendez et al. (2021). These works have allowed for deriving the dimensionless transfer function of the coating film against harmonic and non-harmonic jet oscillations and/or pulsations, showing that a narrow range of high sensitivity exists in the dimensionless frequency domain obtained by scaling the problem using time and length scales built from the liquid film properties. This result is in agreement with the findings of Lunz & Howell (2018) for liquid coatings on inclined planes.

Nevertheless, the main limitation of these theoretical works is that they are based on a one-way coupling formulation in which synthetic perturbations are introduced into the film to analyse its receptivity. More advanced two-way coupling requires more sophisticated simulations and modelling approaches. Examples of such two-way coupling approaches are proposed by Lavalle et al. (2015), who combined an integral film model with a full Navier–Stokes solver for the shearing laminar gas flow; Miao, Hendrickson & Liu (2017), who proposed a similar solution using an immersed boundary (IB) method; and Ojiako et al. (2020), who coupled a single-phase solver with a film model for a gas jet impinging on a liquid bath to obtain the stationary shape of the gas–liquid interface. However, the use of such kind of ‘hybrid’ approaches for the jet-wiping problem remains an open challenge due to the large difference in the scales of the liquid and the (turbulent) gas flow, especially in the wiping conditions encountered in galvanization.

Alternative solutions using state-of-the-art methods from two-phase computational fluid dynamics (CFD) face the same problems. The most common approach is to combine the volume of fluid (VOF) method and large-eddy simulation (LES) modelling for gas jet turbulence. This was the method used by Myrillas et al. (2013) and Barreiro-Villaverde et al. (2021) for the three-dimensional (3-D) simulations reproducing the wiping of highly viscous fluids in the experiments of Gosset (2007), Gosset et al. (2019) and Mendez et al. (2019b), which allowed for the first description of the wiping instability. However, these investigations focused on conditions that are extremely far from the ones encountered in galvanizing conditions, as further discussed in the following section. To the best of the authors’ knowledge, high-fidelity simulations of the process in galvanizing conditions are limited to two-dimensional (2-D) computations (Pfeiler et al. 2017b,a) while the only attempt to reproduce it in three dimensions is due to Aniszewski et al. (2020). Despite the significant amount of computational resources (${\approx }5\times 10^{5}$ central processing unit (CPU) hours) and the use of highly optimized grid-adapting code such as Basilisk (Popinet 2015), the simulations in Aniszewski et al. (2020) only covered a few instants of the initial stage of the interaction, when the jet first impinges on the liquid film and begins to produce the wiping meniscus. Extending those simulations to achieve established wiping conditions was not feasible.

Considering the prohibitive computational costs reported in Aniszewski et al. (2020), the full high-fidelity simulation of the jet wiping in galvanizing conditions is likely going to remain unfeasible in the near future. Therefore, an important open question is the extent to which the results obtained for highly viscous fluids in Barreiro-Villaverde et al. (2021), Gosset et al. (2019) and Mendez et al. (2019b) are relevant for the much broader spectra of wiping conditions encountered in industry. This article is an attempt to tackle this question. Specifically, we complement the work in Barreiro-Villaverde et al. (2021) with the investigation of the wiping process in conditions that are closer to the industrial ones. Although far from achieving full dynamic similarity, as discussed in the following section, the investigated conditions cover a completely different wiping regime than what was previously reported in the literature. Therefore, the results give an insight into the dynamics and the scaling of the jet-wiping instability and perhaps enable an educated extrapolation. Moreover, we use the CFD results to examine the theoretical foundations of the classic jet-wiping models.

The rest of the manuscript is organized as follows. Section 2 presents the investigated conditions and the scaling approaches. Section 3 reports on the methodology, including the numerical approach (§ 3.1), the data-driven modal decomposition (§ 3.2) and the simplified jet-wiping model (§ 3.3) tested against the CFD. Finally, § 4 collects the results and § 5 the conclusions.

2. Selected test cases and scaling laws

We consider planar gas jet wiping, as sketched in figure 1. The slot nozzle of opening $d$ and width $W\gg d$ is located at a distance $Z$ from a vertical strip moving upwards at speed $U_p$. The strip is flat, and the problem is treated as isothermal. The relevant fluid properties are the density $\rho _g$ and kinematic viscosity $\nu _g$ of the gas, the density $\rho _l$ and kinematic viscosity $\nu _l$ of the liquid, and the surface tension $\sigma$ at the interface between the two. The nozzle stagnation chamber is maintained at a gauge pressure $\Delta P_N$, leading to a jet average exit velocity $U_J$. The final coating thickness downstream of the wiping region is $h_f$, and the liquid film forced back by the impingement is usually referred to as run-back flow. The wiping mechanism (that is, the removal of some of the liquid and the consequent thinning of the film) is due to the pressure gradient distribution $\partial _x p_g(x)$ and the shear stress distribution $\tau _g(x)$ produced by the impingement. Both quantities are pictorially illustrated in figure 1 and are often referred to as ‘wiping actuators’.

Figure 1. Schematic of the jet-wiping process, recalling the main operating parameters. The substrate moves upwards at a speed $U_p$, against gravitational acceleration $\boldsymbol {g}$. The typical shape of the streamwise pressure gradient $\partial _x p_g(x)$ and shear stress $\tau _g(x)$ distributions produced by the impinging jet at the film interface are represented on the right-hand side. The average final film thickness is denoted $h_f$.

To analyse the scaling of the problem and normalize the results from the CFD, we introduce a set of reference quantities denoted with square brackets. Therefore, $[a]$ is the reference quantity used to scale the variable $a$, in order to obtain the associated dimensionless variable $\hat {a}=a/[a]$. The chosen set of reference quantities is the one introduced in Mendez et al. (2021) and is listed in table 1 for completeness.

Table 1. Reference quantities used to scale the CFD results in this work.

The film thickness is scaled with respect to $[h] = \sqrt {\nu _l U_p/g}$. This is approximately the maximum thickness that can be withdrawn by the moving strip in the absence of wiping (see Spiers, Subbaraman & Wilkinson 1974). Under the assumption of long wavelengths, the streamwise reference scale is taken as $[x]=[h]/\varepsilon$, with $\varepsilon \ll 1$ a small film parameter. Following Shkadov (1967), we set this parameter to keep the surface tension to the leading order even in the long-wavelength scaling. This yields (see (A6) in Mendez et al. (2021)) $\varepsilon ={Ca}^{-1/3}$ with ${Ca}=\rho _l \nu _l U_p/\sigma$ the capillary number. A natural choice for the reference streamwise velocity $[u]$ is thus the strip velocity $U_p$, since it controls the relative importance of viscosity, inertia and capillary forces in the liquid film. The time scale is accordingly defined as $[t]=[x]/U_p$, while the choice for the remaining reference quantities is in line with the literature on jet wiping (Thornton & Graff 1976; Buchlin 1997) and sets $[p] = \rho _l g [x]$ and $[\tau ] = \mu _l [u]/[h]$.

At a macroscopic level, we follow the scaling discussed by Gosset et al. (2019) and define the wiping conditions in terms of Reynolds numbers $\textit {Re}_f=U_p h_f/\nu _l$ for the liquid film and $\textit {Re}_j=U_j d/\nu _g$ for the gas jet, dimensionless standoff distance $\hat {Z}=Z / d$, Kapitza number $\textit {Ka}=\sigma \rho _l^{-1} \nu _l^{-4/3} g^{-1/3}$ with $g$ the gravitational acceleration and three more numbers linked to the wiping ‘strength’ of the impinging jet in relation to the liquid properties. These are (i) the wiping number ${\varPi }_g= \Delta P_N d/(\rho _l g Z ^2)$, which relates the maximum pressure gradient (responsible for most of the wiping) to the liquid density, (ii) the shear number $\mathcal {T}_g=\Delta P_N d_n / (Z(\rho _l g \mu _l U_p)^{1/2})$, relating the maximum shear stress to the liquid viscosity and (iii) the ‘blockage’ number $[h]/Z$ measuring the intrusiveness of the liquid film on the gas jet flow, to characterize its geometrical confinement.

This article analyses four wiping conditions, for which table 2 collects the aforementioned parameters. Two of these (cases 1 and 2) are taken from Barreiro-Villaverde et al. (2021) and correspond to experiments in Mendez et al. (2019b). These consider the wiping of dipropylene glycol (DG) with a jet of air ($\rho _g=1.2\ {\rm kg}\ {\rm m}^{-3}$ and $\nu _g= 1.48 \times 10^{-5}\ {\rm m}^2\ {\rm s}^{-1}$) with $Z = 18.5$ mm, $d=1.3$ mm and $U_p = 0.34\ {\rm m}\ {\rm s}^{-1}$ and discharge velocity $U_j\approx 26\ {\rm m}\ {\rm s}^{-1}$ (case 1, with $\Delta P_N=425$ Pa) and $U_j\approx 38\ {\rm m}\ {\rm s}^{-1}$ (case 2, with $\Delta P_N=875$ Pa). The other two cases (cases 3 and 4) were added for this work and consider the wiping of water (W) by an air jet with $Z = 10$ mm, $d=1$ mm and $U_p = 1\ {\rm m}\ {\rm s}^{-1}$ and discharge velocity $U_j\approx 42\ {\rm m}\ {\rm s}^{-1}$ (case 4, with $\Delta P_N=1$ kPa) and $U_j\approx 50\ {\rm m}\ {\rm s}^{-1}$ (case 3, with $\Delta P_N=1.5$ kPa). Even if we did not simulate the industrial conditions, it is instructive to compare the selected cases with those on a galvanizing line. Table 2 reports the relevant parameters for a moderate wiping of molten zinc by means of an air jet with $Z = 15$ mm, $d=1.2$ mm and $U_p = 3\ {\rm m}\ {\rm s}^{-1}$, and discharge velocity $U_j\approx 160\ {\rm m}\ {\rm s}^{-1}$ ($\Delta P_N\approx 20$ kPa). These would produce a final thickness of approximately ${\approx }25\ \mathrm {\mu }{\rm m}$.

Table 2. Dimensional and dimensionless wiping conditions for cases 1 and 2 with dipropylene glycol (DG), with $\rho _{l}=1023\ {\rm kg}\ {\rm m}^{-3}$, $\nu _{l}= 7.33\times 10^{-5}\ {\rm m}^2\,{\rm s}^{-1}$, $\sigma _{l} = 0.032\ {\rm N}\ {\rm m}^{-1}$; cases 3 and 4 with water (W), with $\rho _{l}=1000\ {\rm kg}\ {\rm m}^{-3}$, $\nu _{l}= 1\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$, $\sigma _{l} = 0.073\ {\rm N}\ {\rm m}^{-1}$; and an example of galvanizing conditions with zinc (Galvanization), with $\rho _{l}=6500\ {\rm kg}\ {\rm m}^{-3}$, $\nu _{l}= 4.5\times 10^{-7}\ {\rm m}^2\ {\rm s}^{-1}$, $\sigma _{l} = 0.78\ {\rm N}\ {\rm m}^{-1}$.

Table 2 also shows that the newly investigated conditions with water are significantly different from previous test cases and much closer to galvanizing conditions. In particular, the similarity in terms of $[h]/Z$, $h_f/[h]$, ${Ca}$, $\varPi _g$ and $\mathcal {T}_g$ is remarkable, suggesting that the overall wiping mechanism reaches a good degree of similitude. On the other hand, significant differences remain in terms of the Reynolds number in both the liquid and the jet and the Kapitza number, suggesting that the impact of turbulence phenomena in both flows and the interface damping due to surface tension might be underestimated.

3. Methodology

This section describes the numerical configuration of the CFD simulations (§ 3.1), the modal analysis techniques applied to capture the instability mechanism (§ 3.2) and the theoretical foundations of the liquid-film models for jet wiping (§ 3.3).

3.1. Three-dimensional VOF LES simulations

The selected test cases are simulated using the open-source CFD software OpenFOAM (Weller et al. 1998), and the numerical model combines the algebraic VOF formulation of the interFoam solver with the Smagorinski LES model for the turbulence treatment of the gas jet flow. Barreiro-Villaverde et al. (2021) previously validated this approach with experiments, both in terms of stationary (mean final coating thickness) and non-stationary flow features (coherent structures and frequency content of both gas jet and liquid film). The LES technique involves the low-pass filtering of the equations with a Gaussian spatial filter $\boldsymbol {G}_\varDelta (\boldsymbol {x})$ such that $\bar {\phi } = \boldsymbol {G}_\varDelta (\boldsymbol {x}) \circledast \phi$ where $\phi$ denotes the flow variables, the $\circledast$ the convolution operator and the overbar refers to the filtered flow variables. Therefore, the incompressible Navier–Stokes system reads

(3.1a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\bar{{\boldsymbol{u}}} = 0 , \end{gather}$$

(3.1b)$$\begin{gather}\frac{\partial \left( \bar{\rho} \bar{\boldsymbol{u}}\right) }{\partial t} + \left( \bar{\rho} \bar{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla}\right) \bar{\boldsymbol{u}} ={-} \boldsymbol{\nabla} \bar{p} + \bar{\rho} \boldsymbol{g} + \boldsymbol{\nabla} \boldsymbol{\cdot} [ 2\bar{\mu} \bar{\boldsymbol{\mathsf{S}}}] + \bar{F}_\sigma - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} , \end{gather}$$

where ${\boldsymbol {u}}=(u,v,w)$ is the velocity, $p$ the pressure, $\boldsymbol {g}$ is the gravitational acceleration, $\mu$ the dynamic viscosity, $\boldsymbol{\mathsf{S}}=(\boldsymbol {\nabla }{\boldsymbol {u}} + \nabla ^{T}{\boldsymbol {u}})/2$ the strain rate tensor, $F_\sigma$ the term accounting for the Laplace pressure computed using the continuum surface force model of Brackbill, Kothe & Zemach (1992) and $\boldsymbol {\tau }$ the subgrid stress term. This last term accounts for the turbulent stresses at scales smaller than the spatial filter, which in the Smagorinski model are computed as

(3.2a,b)\begin{equation} {\tau_{u} ={-}2 \mu_t \bar{\boldsymbol{\mathsf{S}}} ,} \quad {\mu_t = \rho (C_s \Delta x_{cell})^2 |\bar{\boldsymbol{\mathsf{S}}}| .} \end{equation}

The system is completed with the algebraic VOF implementation in interFoam (Deshpande, Anumolu & Trujillo 2012), which solves an additional transport equation for the liquid volume fraction $\alpha$,

(3.3)\begin{equation} \frac{\partial \bar{\alpha}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\bar{{\boldsymbol{u}}} \bar{\alpha}) + \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\bar{{\boldsymbol{u}}}_r \bar{\alpha} (1-\bar{\alpha}) \right) = 0 , \end{equation}

where the last term introduces an artificial velocity $\bar {{\boldsymbol {u}}}_r$ at the interface to mitigate its smearing. It is important to note that we neglect the extra subgrid terms due to unresolved surface tension terms and surface deformation due to the phase-averaging of VOF and the low-pass filtering of the LES (see Lakehal (2018) for more details). The large-eddy interface simulation models aim at providing closure models for the interfacial subgrid terms, although these are very case dependent and their implementation is only justified in cases where the subgrid scales are relevant to the flow, such as in jet atomization.

The domain and the mesh are shown in figure 2. The domain includes the coating bath from which the flat substrate is withdrawn, and a portion of the slot nozzle from which the gas jet is released. It spans 50 mm above and below the jet axis in the streamwise direction ($x$), 20 mm in addition to the standoff distance $Z$ in the cross-stream direction ($y$) and 20 mm in the spanwise direction ($z$). The boundary conditions are indicated in figure 2, and the gas jet is established through an inlet condition with prescribed stagnation pressure $\Delta P_N$. Turbulence is promoted by introducing white noise perturbation with an amplitude of 5 % of the jet inlet velocity. The lateral patches (not shown) are set to cyclic.

Figure 2. Numerical domain, boundary conditions and mesh discretization.

The discretization is particularly challenging because of the multiscale nature of the problem, requiring a resolution of the order of the micrometre in different regions of the domain and along different directions. Since a complete mesh sensitivity analysis was out of the reach of our computational resources, the mesh design was based on (i) scaling considerations, (ii) the successful experimental validation of the simulations in Barreiro-Villaverde et al. (2021) and (iii) a set of preliminary tests. More specifically, the meshing near the moving strip was carried out keeping the same aspect ratio and the same dimensionless resolution across the liquid film as in Barreiro-Villaverde et al. (2021).

These considerations lead to a cell size in the region of gas–liquid interaction ($-Z< x< Z$) of $\varDelta _x=50\ \mathrm {\mu }{\rm m}$ and $\varDelta _y=2\ \mathrm {\mu }{\rm m}$ across the film thickness. This provides 10 to 15 mesh points within the liquid film and 100 to 200 mesh points within one wavelength of the leading mode in the streamwise direction. This wavelength was first estimated from the results in Barreiro-Villaverde et al. (2021), scaled according to the reference quantities in table 1, and then confirmed a posteriori by the modal decomposition of the results. In dimensionless coordinates, this is much finer than the mesh in Aniszewski et al. (2021) which allowed for a maximum of four points within the liquid film in the 2-D simulation and less than one mesh point within the film in the 3-D simulations.

On the gas side, the mesh was designed using the Taylor microscale in the nozzle as a reference. This characteristic length was estimated as $\lambda _T = (10 k \nu _g / \epsilon )^{0.5}$, with $k$ the turbulent kinetic energy and $\epsilon$ the turbulent dissipation, computed using the empirical procedure in Rodriguez (2019),

(3.4a,b)\begin{equation} {k=\frac{3}{2} (U_J I )^2, \qquad \epsilon =C_{\mu} \frac{k^{2/3}}{ \ell },} \end{equation}

with $U_J$ the mean velocity of the jet, $I$ the turbulence intensity estimated as ${I=0.144 Re_j^{-0.146}}$, $C_\mu = 0.09$ and $\ell =0.07 D_h$ the integral length scale based on the hydraulic diameter $D_h \approx 2d$. This estimation leads to $\lambda _T\approx 230\ \mathrm {\mu }{\rm m}$ for case 3 and $\lambda _T\approx 220\ \mathrm {\mu }{\rm m}$ for case 4. The mesh resolution in the streamwise direction was thus set to $\varDelta _y = 100\ \mathrm {\mu }{\rm m}$, as a good compromise between resolution and cost and to keep $\varDelta _y/\lambda _T<1/2$ in all simulations. In absolute terms, the resolution is comparable to the one in Aniszewski et al. (2021). However, these authors consider a flow at much higher velocity (${{Re}_J\approx 14\,300}$), leading to $\varDelta _y/\lambda _T>1$. In the cross-streamwise direction, the resolution was set to $\varDelta _x = 10\ \mathrm {\mu }{\rm m}$, leading to $y^+<5$ at the nozzle's wall for all cases. Finally, the spanwise resolution was set to $\varDelta _z=250\ \mathrm {\mu }{\rm m}$, keeping a similar aspect ratio in the spanwise direction as in Barreiro-Villaverde et al. (2021) and expecting a predominantly 2-D dynamics in the flow as in the previous works. The final mesh count reaches $14\times 10^6$ hexahedral cells.

The mesh quality was also assessed a posteriori using the LES index proposed in Celik, Cehreli & Yavuz (2005),

(3.5)\begin{equation} L_e = \frac{1}{1+\alpha_{\nu} \left( \dfrac{\nu_t}{\nu} \right)^n} , \end{equation}

with $\alpha _{\nu } = 0.05$ and $n=0.53$ from correlations, and $\nu _t=\mu _t/\rho _g$ the turbulent viscosity from the simulations, computed from (3.2a,b). The LES index is above $80\,\%$ in the wiping region in all cases, following the recommendations in Pope (2000).

The simulations are initialized as in Barreiro-Villaverde et al. (2021) following two steps. The first phase begins with a uniform thickness distribution on the strip with $\hat {h}=1$ and ends with the establishment of a drag-out flow in the absence of the jet. This requires approximately $100$ ms of physical time. The second phase introduces the jet flow and ends after the wiping meniscus is established. This requires approximately 100 ms to simulate the initial transient, after which 350 ms of physical time is computed. The time step is $10^{-7}$ s to keep the Courant–Friedrichs–Lewy number below 0.9, and a total of $n_t=3500$ instantaneous flow fields are exported with a sampling rate of 10 kHz. The computation requires approximately 1000 hr running in parallel on 512 Intel E5–2680v3 CPUs from the Centro de Supercomputacion de Galicia (CESGA) for one test case with water (case 3), and a total of $5\times 10^6$ CPU hours for the four test cases.

3.2. Multiscale proper orthogonal decomposition

The results are processed using the multiscale proper orthogonal decomposition (mPOD) (Mendez, Balabane & Buchlin 2019a; Mendez 2023) and its extension, as in Barreiro-Villaverde et al. (2021). This technique decomposes the data $\boldsymbol{\mathsf{D}}(\boldsymbol {x},t)$ as a linear combination of modes such that

(3.6)\begin{equation} {\boldsymbol{\mathsf{D}}(\boldsymbol{x},t)=\sum^{R}_{r=1} \sigma_r \boldsymbol{\phi}_r(\boldsymbol{x}) \boldsymbol{\psi}_r(t),} \end{equation}

where $\sigma _r$ is the mode amplitude, $\boldsymbol{\phi} _r(\boldsymbol {x})$ is the mode spatial distribution and $\boldsymbol{\psi}_r (t)$ its temporal evolution. The structures $\boldsymbol{\phi} _r$ and $\boldsymbol{\psi} _r$ form, respectively, a basis for spatial distribution and temporal evolution of the data, while $\sigma _r$ indicates the energetic weight of each structure in the dataset. The mPOD is a data-driven decomposition that bridges the gap between the energy-based proper orthogonal decomposition (POD) (Sirovich 1991; Holmes, Lumley & Berkooz 1996) and frequency-based dynamic mode decomposition (Rowley et al. 2009; Schmid 2010) methods. In particular, the mPOD modes are optimal modes for a prescribed frequency partition and have a band-limited frequency content.

In a first step, we decompose the normalized film thickness, defined as $\check {h} (x,z,t) = ({h(x,z,t) - \bar {h}(x,z)})/{\sigma _h (x,z)}$ with $\bar {h} (x,z)$ the average thickness and $\sigma _h (x,z)$ the thickness standard deviation (Barreiro-Villaverde et al. 2021). This normalization allows the decomposition to equalize the importance of waves upstream and downstream of the wiping region despite the largely different thickness. Since we are interested in the most dominant travelling wave patterns, described by two sinusoidals in quadrature (phase-shifted by ${\rm \pi} /2$) both in time and space, we only retain the pair of modes $\mathcal {M}$ with the largest and equal amplitudes ($\sigma _{\mathcal {M}}$), having spatial ($\boldsymbol{\phi}_{\mathcal {M}}$) and temporal structures ($\boldsymbol{\psi}_{\mathcal {M}}$) in quadrature.

In a second step, we extract the most correlated coherent structures in the gas jet by projecting variables linked to the jet flow (e.g. velocity field, pressure or shear stress distribution at the wall) onto the dynamics of the leading wave patterns. Therefore, denoting as $\boldsymbol{\mathsf{D}}_g (\boldsymbol {x},t)$ the general jet flow variable and as $\boldsymbol{\psi} _{\mathcal {M}}(t)$ the temporal structures of the leading wave pattern in the film thickness, the spatial distribution of the most correlated pattern is defined as

(3.7) \begin{equation} {\boldsymbol{\phi}_g(\boldsymbol{x}) = \frac{\langle \boldsymbol{\mathsf{D}}_g (\boldsymbol{x},t), \boldsymbol{\psi}_{\mathcal{M}} (t)\rangle_t}{\sigma^2_{\mathcal{E}}}\quad \mbox{with}\ \sigma_{\mathcal{E}}=|| \langle \boldsymbol{\mathsf{D}}_g (\boldsymbol{x},t), \boldsymbol{\psi}_{\mathcal{M}} (t)\rangle_t||_s} , \end{equation}

where $\langle \boldsymbol{a}(t),\boldsymbol{b}(t)\rangle _t=1/T\int ^T_0 \boldsymbol{a}(t) \boldsymbol{b}(t)\,{\rm d}t$ is the inner product in time, and $||\boldsymbol{a}(\boldsymbol {x})||^2_s=1/|\varOmega |\int _{\varOmega } \boldsymbol{a}^2(\boldsymbol {x})$ is the $L_2$ norm in space, with $\varOmega$ the simulation volume and $|\varOmega |$ its measure. The denominator in (3.7) has the scope of having $||\boldsymbol {\phi} _{g}(\boldsymbol {x})||_s=1$. The projection of a dataset (e.g. velocity field) onto the leading structures of another dataset (e.g. film height evolution) was first proposed by Borée (2003) for the POD and referred to as ‘extended’ POD. Barreiro-Villaverde et al. (2021) use the same idea for the ‘extended’ mPOD (emPOD).

3.3. One-dimensional analytical model for jet wiping

Analytical models for jet wiping were initially derived to predict the final thickness profiles as a function of the process parameters (Thornton & Graff 1976; Tuck 1983; Ellen & Tu 1984; Yoneda 1993; Buchlin 1997), and then extended to quasisteady formulations by Hocking et al. (2011) and more advanced integral formulations in Mendez et al. (2021).

Early validations of these models were limited to the final film thickness. A more extensive validation of the whole average thickness distribution for steady state models was presented by Lacanette et al. (2006), while Barreiro-Villaverde et al. (2021) explored the validity of integral models in dynamic conditions. However, both Lacanette et al. (2006) and Barreiro-Villaverde et al. (2021) considered the case of highly viscous fluids. In this work, we extend the validation to the wiping of water, both in time averaged and time-dependent conditions.

Thin film models build upon the 2-D incompressible Navier–Stokes equations for the liquid film, simplified using the longwave approximation. The later assumes that the characteristic length scale in the streamwise direction ($x$) is much larger than in the cross-stream ($y$), i.e. the perturbation wavelength $\lambda$ is large compared with the liquid film thickness $h$. The separation of scales is provided by the film parameter $\varepsilon = [h]/[x]=Ca^{1/3}$ introduced in § 2. Neglecting terms of order $>O(\varepsilon )$, integrating the equations along the cross-stream direction using the Leibniz rule yields the following system of PDEs constituting the IBL model (Mendez et al. 2021):

(3.8a)$$\begin{gather} {\partial_{\hat{t}} \hat{h}+\partial_{\hat{x}} \hat{q}=0, } \end{gather}$$

(3.8b)$$\begin{gather}{\varepsilon \,Re \left( \partial_{\hat{t}} \hat{q}+\partial_{\hat{x}} \mathcal{F} \right)=\hat{h} \left( 1-\partial_{\hat{x}} \hat{p}_{g}+\partial_{\hat{x} \hat{x} \hat{x}} \hat{h}\right) + \hat{\tau}_{g} - \hat{\tau}_{w} ,} \end{gather}$$

where ${Re}=[u][x]/\nu _l = \sqrt {U_p^3 /g \nu _l}$ is the Reynolds number, $\hat {q}$ is the volumetric flow rate per unit width, $\partial _{\hat {x}} \hat {p}_g$ and $\hat {\tau }_g$ are the pressure gradient and shear stress distribution at the interface produced by the impinging jet, $\mathcal {F} = \int _{0}^{\hat {h}} \hat {u}^{2} \,{\rm d} \hat {y}$ is the integral advection term and $\hat {\tau }_{w} = \partial _{\hat {y}} \hat {u}|_{\hat {y}=0}$ is the shear stress at the wall. These last two terms require some assumptions on the streamwise velocity profile. The simplest IBL model assumes a parabolic, self-similar profile. Using the boundary conditions and the definition of the flow rate per unit width, it can be written as

(3.9)\begin{equation} \hat{u}= \left(\frac{3 \hat{q}}{\hat{h}}-\frac{3}{2} \hat{h}+3 \right) \left[\left(\frac{\hat{y}}{\hat{h}}\right)-\frac{1}{2}\left(\frac{\hat{y}}{\hat{h}}\right)^{2}\right]+\hat{\tau}_{g} \hat{y}-1, \end{equation}

leading to the following closure for (3.8):

(3.10a,b)\begin{equation} {\mathcal{F}=\frac{\hat{h}^{3} \hat{\tau}_{g}^{2}}{120}+\frac{\hat{h} \hat{q} \hat{\tau}_{g}}{20}+\frac{\hat{h}^{2} \hat{\tau}_{g}}{20}+\frac{6 \hat{q}^{2}}{5 \hat{h}}+\frac{2 \hat{q}}{5}+\frac{\hat{h}}{5}} \quad{\text{and}}\quad \tau_w ={-}\frac{3\hat{q}}{\hat{h}^2} - \frac{3}{\hat{h}}. \end{equation}

It is worth noticing that the boundary condition is $\hat {u}(0)=-1$: keeping the classic convention in the literature of falling liquid films, the streamwise axis points downwards (same direction as gravity), contrary to the axes in the numerical simulations presented in this work.

The integral model in Hocking et al. (2011) can be obtained by neglecting inertia, hence $\mathcal {F}=0$. In a further simplification, the steady state models in Tu & Ellen (1986) and Buchlin (1997) can be reached by also removing the surface tension and time dependent terms. In these conditions, since $\partial _{\hat {x}} \hat {q} = 0$, the system of equations in (3.8) reduces to a cubic equation in terms of film thickness $\hat {h}$, with $\hat {q}$ left as a free parameter,

(3.11)\begin{equation} {\left( \partial_{\hat{x}} \hat{p} + 1 \right) \hat{h}^3 - \tfrac{3}{2}\hat{\tau}_g \hat{h}^2 - 3\hat{h} + 2 \hat{q} = 0}. \end{equation}

This is known as the zero-order model. At each position $x$, the cubic equation (3.11) admits two positive solutions. The blending between these allows recovering the full interface distribution. Independently from their respective degree of approximation, all these models require as inputs the distributions of pressure gradient and shear stress at the film interface (see Gosset (2007) for a review). This is why to assess these models in the best possible conditions, we extract the stress distributions from the numerical results.

4. Results

The results of the two-phase simulations are first compared with thin film models in terms of average thickness and velocity profiles in § 4.1. Section 4.2 analyses the wave patterns on the liquid film, while § 4.3 focuses on the gas flow structures that are most correlated with the leading wave patterns. Finally, § 4.4 discusses the mechanism of interaction between the two phases. In all sections, the new results on cases 3 and 4 are analysed along with the ones previously obtained for cases 1 and 2. The spatiotemporal distribution of thickness $h(x,z,t)$ and streamwise pressure gradient $\partial _x p(x,z,t)$ of cases 1 and 3 are openly available at http://doi.org/10.17605/OSF.IO/QAGBZ.

4.1. Cross-assessment of CFD and liquid-film models

As explained previously, the distributions of pressure gradient and shear stress at the interface are required to assess the thin film models in (3.3). These were retrieved from the postprocessing of the numerical results.

Starting from the averaged predictions, figure 3(a) shows the time-averaged distribution of pressure at the interface in the middle plane, i.e. $z=L_z/2$, for case 3. A Gaussian distribution $\hat {p}_g(x)=\hat {p}_M\exp (-\hat {x}^2/\gamma _p)$ with $\hat {p}_M=43.5$ and $\gamma _p = 0.992$ appears appropriate. Figure 3(b) shows the link between the time averaged pressure gradient and the time averaged shear stress at the interface. Within the range $\hat {x}\in [-1,1]$, a reasonable correlation is provided by the linear trend $\tau _g(x)\approx -0.07 \partial _{\hat {x}}\hat {p}(x)+0.5\,\mbox {sign}(x)$, with $\mbox {sign}(x)=1$ if $x>0$, $\mbox {sign}(x)=-1$ if $x<0$ and $\mbox {sign}(0)=0$. The same approach was used for case 4.

Figure 3. Time-averaged dimensionless pressure distribution at the interface and the result of exponential fitting (a), and correlation between the pressure gradient distribution and the shear stress distribution at the interface (b). In both plots, the blue continuous line represents the CFD data, while the dashed red line represents the Gaussian and the linear fitting, respectively.

Introducing these distributions in the models (3.8) and (3.11) yields the results in figure 4. Figure 4(a) compares the results in terms of mean thickness profile for the CFD (orange solid line), the full IBL model in steady conditions (red dashed line) and the zero-order model (blue dotted line). The figure is complemented with a close-up above the wiping region and shaded areas to show the envelope of film thickness distributions. These were filtered using the emPOD to focus only on the leading travelling patterns. The model predictions were obtained using the pressure gradient and the shear stress distributions at the interface obtained by fitting an exponential curve on the pressure distribution and a linear fit between this and the shear stress. These are compared with the envelope values after filtering these quantities with the emPOD in figures 4(b) and 4(c), respectively.

Figure 4. Comparison between the CFD time averaged thickness profiles at $z=L_z/2$ and the theoretical predictions using the zero-order and IBL model (a), together with the envelopes of the pressure gradient $\partial _x p$ (b) and shear stress $\tau _x$ (c). The solid lines in figures 4(b) and 4(c) represent the inputs of the theoretical models and have been obtained via Gaussian regression.

The overall agreement between the thin film models and the CFD is satisfactory in terms of average film thickness, especially considering the simplifying assumptions underpinning the models and the uncertainty on the wiping actuators. The IBL model provides an excellent prediction of the final film thickness and a satisfactory prediction of the meniscus shape. The discrepancy between the IBL and the zero-order model suggests that the surface tension and the advection term $\mathcal {F}$, both neglected in the latter, are important. The large fluctuations of both the film thickness and wiping actuators illustrate the importance of the flow unsteadiness, analysed in detail in §§ 4.3 and 4.4.

Moving now to the dynamic conditions, the main interest is to analyse the validity of the self-similar parabolic profile underpinning the IBL model. This was shown to be surprisingly valid for cases 1 and 2 in Barreiro-Villaverde et al. (2021), despite the value of the film parameter ($\varepsilon =[h]/[x]=0.93$), certainly too close to unity to justify the longwave approximation. Moreover, these cases are characterized by ${Re}\approx 7.4$, leading to a reduced Reynolds $\delta \approx \varepsilon \,{Re}= 6.8$ which is also outside the theoretical limit $\delta \sim O(1)$ above which the integral boundary layer approximation is in principle valid, as detailed in Kalliadasis et al. (2012). Cases 3 and 4 are characterized by a lower $\varepsilon =0.24$ but a much higher ${Re}\approx 320$, leading to $\delta \approx 77$. Figure 5 shows for case 1 (figure 5a) and case 3 (figure 5b), a contour plot of the streamwise velocity component in the $z$-midplane ($u(x,y,z=L_z/2)$) on the left-hand, with the horizontal red lines indicating the location from which the velocity profiles are extracted and plotted on the right-hand. These plots compare the CFD data (dashed blue) with the self-similar profile in (3.9). All variables are scaled by the reference quantities in table 1.

Figure 5. Contour map of the streamwise velocity with streamlines at the $z$-midplane. The data are shown both in dimensional (left-hand side) and dimensionless axes (right-hand side). The horizontal red lines in the contours indicate the locations of the velocity profiles on the right-hand side. The continuous lines show the profile obtained with the simulation, and the dashed lines the self-similar profile in (3.9). (a) Case 1; (b) case 3.

The comparison of five profiles at different positions in the streamwise direction shows that the velocity profile remains close to parabolic, even at the much higher Reynolds number in case 3. This is true in all the domain for case 1 and only in absence of 3-D effects in case 3. For the illustrated snapshots, the validity of the parabolic assumption is particularly striking in the profile at $x \approx -4$ mm ($\hat {x}<-2.2$) for case 1 and the profile at $x \approx -1.8$ mm ($\hat {x}=-1.3$) for case 3. Both are characterized by a large interface curvature and shear stress (see also figure 4c). Within $\hat {x}\in [-1,1]$, the dynamics of the flow field within the liquid film is quite similar, while differences become more pronounced as the waves in the run-back flow develop. At $x<-2$ mm ($\hat {x}<-1.7)$), the flow becomes 3-D for case 3 (as further illustrated in the following) and vortices appear between the crests of two successive waves. Therefore, while the parabolic assumption is appropriate in a large portion of the domain for case 1, its validity is confined to $\hat {x}>-1.5$ for case 3. Nevertheless, as discussed in the following sections, this is the area where most of the interaction between the gas jet and the liquid film takes place. Finally, we assess the validity of the key ingredient of the boundary layer approximation, allowing for neglecting the convective/inertial effect in the cross-stream direction in the film, and assuming zero cross-stream pressure gradient. To this end, figure 6 shows the pressure profiles at the interface $\hat {p}_i=\hat {p}(\hat {y}=\hat {h})$ and at the wall $\hat {p}_w = \hat {p}(\hat {y}=0)$ (figure 6a), the instantaneous film thickness distribution (figure 6b) and the standard deviation of the differential pressure across the film $\sigma _{\Delta \hat {p}}$, where $\Delta \hat {p} = \hat {p}_i - \hat {p}_w$ (figure 6c). The pressure at the wall and at the interface are identical in the final film, while some variations are visible in the run-back flow, where the flow is forced to reverse direction due to the wiping effect of the jet. This, together with the larger thickness and local Reynolds number make the cross-stream inertia more important and challenge the validity of the boundary layer approximation.

Figure 6. Dimensionless pressure profiles at the wall $\hat {p}_w = \hat {p}(\hat {y}=0)$ and at the interface $\hat {p}_i = \hat {p}(\hat {y}=\hat {h})$ (a) and thickness (b) for a given time step, together with the standard deviation $\sigma _{\Delta \hat {p}}$ of the dimensionless pressure difference across the film $\Delta \hat {p} = \hat {p}_w - \hat {p}_i$ (c).

4.2. Wave patterns on the liquid film

Figure 7 illustrates the dynamics of the liquid film upstream of the wiping region ($x<0$) in case 1 (figure 7a) and 3 (figure 7b). In each plot, an instantaneous film thickness distribution is shown, with the liquid going upwards in red ($u>0$) and the liquid going downwards in green ($u<0$). In case 1, the liquid film is clearly 2-D and features large waves originating approximately $2.5$ mm below the jet axis at a frequency of approximately $20$ Hz, corresponding to $\hat {f}=f{[t]} \approx 0.11$ in the scaling from table 1. These falling waves have a wavelength of the order of $\lambda \approx 12$ mm and evolve over an average film thickness of the order of $h_r\approx 2h_0=3$ mm. The local Reynolds number in this region, defined as ${Re}_r=h_r \Delta U/\nu$, with $\Delta U = \max (U_p - u|_{y=h}) \approx 2U_p$ is of the order of $Re_r\approx 18$. The flow is fairly laminar, and the 2-D nature of the falling waves is in line with what happens at the onset of interface instability in falling liquid films. Case 3 is clearly in a different wiping regime, especially from the liquid perspective as indicated in table 2. The frequency of wave formation ($f \approx 100\unicode{x2013}150$ Hz) in the introduced scaling is $\hat {f} \approx 0.15\unicode{x2013}0.2$, a value notably close to the one in case 1. The reasonable matching of the dimensionless frequencies in the two wiping regimes reinforces the choice of the reference quantities used for scaling, and justifies the dimensionless representation adopted in subsequent comparisons. The effect of increasing the film Reynolds number in case 3 ($Re_r\approx 800$) is clearly visible on the wave patterns on the run-back flow. The waves emerge 2-D – similar to case 1 – at approximately $0.5$ mm below the jet axis, but they rapidly undergo a transition towards a 3-D pattern, reminiscent of what occurs in falling liquid films (see Demekhin et al. 2007).

Figure 7. Snapshots of the film thickness dynamics in cases (a) 1 and (b) 3. The liquid is coloured in red where $u>0$ and in green where $u<0$. A supplementary movie of the 3-D reconstruction of the liquid film is also provided for both cases and is available at https://doi.org/10.1017/jfm.2024.553.

Moving to the multiscale modal analysis of the liquid film thickness, the convergence of the decomposition is illustrated in figure 8 as a function of the normalized energy content of the modes $(\sigma _r / \sigma _1)^2$. The modes denoted with a full marker correspond to the ones that are selected as dominant wave patterns following the criterion defined in § 3.2. For each test case, figure 9 shows a snapshot of the normalized liquid film thickness (figure 8a i,b i,c i,d i), its reconstruction using the leading mPOD modes linked to the dominant travelling wave pattern (figure 8a ii,b ii,c ii,d ii) and the spectra of the associated temporal structures (figure 8a iii,b iii,c iii,d iii). The wiping conditions are recalled in the caption. The spectra computed from the temporal structure of the leading modes are represented as a function of both the dimensional $f$ and dimensionless $\hat {f}$ wave frequencies. Note that the stream and spanwise coordinates are shown in dimensionless form to facilitate the comparison between cases, since the dimensional wavelengths (and frequencies) are very different, as illustrated previously with figure 7.

Figure 8. Convergence of the mPOD decomposition for all cases. The dominant modes in figure 9 are denoted with a full marker.

Figure 9. Comparison between snapshots of the normalized film thickness $\check {h}(\hat {x},\hat {z})$ (i) with the leading wave patterns detected via mPOD (ii). The spectra of the temporal structure of these modes are shown in (iii) as a function of the dimensional $f$ and dimensionless $\hat {f}$ wave frequencies. (a) Case 1 (DG), $\hat {Z}=14.2$, ${\varPi _g =0.16}$, $\mathcal {T}_g = 0.24$; (b) case 2 (DG), $\hat {Z}=14.2$, $\varPi _g =0.33$, $\mathcal {T}_g = 0.41$; (c) case 3 (W), $\hat {Z}=10$, $\varPi _g =1.02$, ${\mathcal {T}_g = 2.23}$; (d) case 4 (W), $\hat {Z}=10$, $\varPi _g =1.53$, $\mathcal {T}_g = 3.02$.

The slow convergence of the decomposition illustrates the complexity of the flow, as it appears from the normalized film thickness $\check {h}$ distributions on the left-hand side. Even if the original (unfiltered) fields depict different dynamics, the mPOD analysis detects similar dominant wave patterns in both the run-back flow and final film. In all the cases, the leading wave pattern is remarkably 2-D, despite the vastly different wiping conditions, and regardless of the 3-D patterns arising upstream of the wiping region with water.

This is linked to the fact that the energy of the modes is more biased towards the final film in the cases with water because the 2-D dominant pattern, although present, is less important in the run-back flow. The energy in the run-back flow is distributed within more modes, featuring 3-D waves at slightly different wavelengths and frequencies. The dominant wavelength upstream and downstream of the wiping differs because of the different advection velocities in the two regions, but like in the run-back film, the range of dimensionless frequency is surprisingly similar (0.1–0.3). In line with the experimental findings with low $\textit {Ka}$ liquids in Gosset et al. (2019), the frequency tends to increase with the wiping strength, and the gas jet Strouhal number ${St}=f Z / U_j$ based on the dominant frequency of the waves, is in the range (0.01–0.04), well below the typical values in free-jet instability and hydrodynamic feedback mechanisms (Mendez, Scelzo & Buchlin 2018). These results highlight the important coupling between the gas jet and the liquid. The following section analyses the gas jet structures linked to those dominant wave patterns.

4.3. Gas jet structures correlated with the leading waves in the liquid film

The dynamics of the gas jet and its connection to the dominant wave patterns has seldom been explored in the literature, due to the challenges in the experimental characterization and the computational cost of high-fidelity simulations. Figure 10 depicts the complexity of the gas jet dynamics for case 1 (figure 10a) and 3 (figure 10b). For both cases, the velocity contours at the $z$-midplane are shown on the left-hand side while a 2-D view of the film in a plane parallel to the moving substrate is shown on the right-hand side. The gas jet structures, extracted using the $Q$-criterion, are shown in grey while the liquid film interface is shown in represented in blue shades.

Figure 10. Flow visualization of case (a) 1 and (b) 3. The gas jet vortices extracted using the $Q$-criterion (grey) are superimposed to the liquid film interface (blue) and a contour plot of the velocity magnitude in the (a) $x\unicode{x2013}y$ plane and (b) $x\unicode{x2013}z$ plane for each test case: (a) case 1 (DG); (b) case 3 (W).

The dynamics of the jet are comparable, as expected from the similarity of the jet Reynolds number ${Re}_j$ (see table 2). At a short distance from the nozzle, 2-D vortices are generated in the shear layers, due to the velocity gradient between the jet core and the quiescent atmosphere. At impingement, these vortices produce a low-amplitude/high-frequency flapping of the jet tip. Then, the main flow splits into two ‘side jets’ that remain attached to the gas–liquid interface until the run-back waves become too steep and trigger flow detachment. Downstream of the wiping, the 2-D vortices are convected along with the side jets and in close interaction with the interface. Finally, these structures undergo transverse modulation and dissipate into smaller eddies. Upstream of the wiping, it is interesting to note that the 2-D–3-D transition seems to be enhanced by the presence of the waves on the run-back flow, especially in case 3 where the vortices become 3-D very close to the impact. Comparatively, they remain unaltered for longer distances on the final film, which suggests the gas–liquid interaction is negligible in this region. These flow visualizations confirm that the gas jet perturbation is essentially 2-D, and therefore, uncorrelated with the formation of 3-D waves on the liquid, as expected from the analysis in § 4.2. This was previously observed in Barreiro-Villaverde et al. (2021) for case 2 using correlation analysis, and is confirmed here for different wiping conditions.

We now focus on the spatial structures in the gas jet flow using the emPOD, to reveal the flow structures that are most correlated with the leading wave patterns on the liquid film (described in figure 9). Figure 11 shows two representative snapshots of the gas velocity field $\boldsymbol {u}=(u,v)$ at the $z$-midplane and its projection for case 1 (figure 11a) and case 3 (figure 11b). Cases 2 and 4 lead to similar results and are thus omitted for brevity, but an animation of the film thickness evolution and gas jet flow for each test case is provided in the supplementary material. In each figure, the contour on the left-hand side is the mean centred flow field (i.e. $\boldsymbol {u}^\prime = \boldsymbol {u}(\boldsymbol {x}) - \bar {\boldsymbol {u}}$), and the one on the right-hand side is the projection of the flow field on the leading mPOD modes of the normalized film thickness.

Figure 11. Mean-centred $u^\prime$ and emPOD projected fields for cases (a) 1 and (b) 3. The plots are scaled such that the horizontal axis ranges between $0< y/Z<1$, and the vertical axis between $-10< x/d<10$. A supplementary movie for each test case in table 2 is also provided. (a) Case 1 (DG); (b) case 3 (W).

While the original flow fields reveal the full spectrum of flow structures, the emPOD filtered fields retain only those that are temporally correlated with the main patterns on the liquid layer or, from a physical perspective, those that are mostly correlated to the gas–liquid coupling. For instance, the vortices in the shear layers (label ‘$v_s$’) and their downstream evolution (label ‘$v_f$’) are not visible in the emPOD counterpart, confirming that these structures are not correlated with the main wave pattern. Their impact on the liquid film is limited to small-scale disturbances that are rapidly damped due to viscosity and capillarity (Barreiro-Villaverde et al. 2023). Here, the leading structures consist of a larger vortex induced by the falling liquid waves (label ‘$v_r$’) and a rigid oscillation of the impinging jet (label ‘$d$’). In the case of low Kapitza liquids, this mechanism was extensively documented in Myrillas et al. (2013), Mendez et al. (2019b) and Barreiro-Villaverde et al. (2021), where a hypothesis on their interaction with the liquid film dynamics was formulated. In particular, it was postulated that the jet oscillation was linked to the vortex–liquid film interaction. The current results with high $\textit {Ka}$ liquids show for the first time that this interaction is much weaker in these conditions because of the comparatively smaller thickness of the liquid film; yet, these structures ‘$v_r$’ and ‘$d$’ persist. Is it because the gas–liquid interaction is similar despite the different flow conditions? Which comes first, the waves in the liquid or the oscillations/deflections in the gas? Is it possible to extrapolate these results to industrial galvanization? These questions are addressed in § 4.4.

4.4. Dynamics of the gas–liquid instability

Previous sections focused on the relevant flow patterns in the liquid film and the gas jet; this section investigates their impact on the momentum transfer at the interface. For two snapshots separated by half-undulation period $T$, figure 12 shows the emPOD filtered dimensionless pressure gradient $\partial _{\hat {x}} \hat {p} = \partial _x p [x]/[p]$ and dimensionless shear stress $\hat {\tau }_{\hat {x}} = \tau _x / [\tau ]$ for case 1 (figure 12a) and case 3 (figure 12b).

Figure 12. Spatiotemporal evolution of the pressure gradient $\partial _{\hat {x}} \hat {p}$ and shear stress $\hat {\tau }_{\hat {x}}$ in case (a) 1 and case (b) 3 for two time steps separated by half an undulation cycle $T$. (a) Case 1 (DG); (b) case 3 (W).

In both cases, the distributions downstream of the wiping ($\hat {x}>0$) remain more or less unaffected – as expected from the negligible activity shown in figure 11 at these locations. The unsteadiness of the wiping actuators is concentrated in the run-back flow, featuring pulsations and oscillations of the wiping point $x^*$ ($\partial _x p (x^*) = \max (\partial _x p)$) and the shear point $x_\tau$ ($\tau _x (x_{\tau }) = \max ({\rm abs}(\tau _x)$), due to the interaction between the lower side jet and the run-back waves.

The connection between the dynamics of the liquid film and wiping actuators is analysed by means of the spatiotemporal contour of the normalized film thickness $\check {h}$ in figure 13 for case 1 (figure 13a) and 3 (figure 13b). This contour map is complemented with the temporal evolution of (i) the impingement point $x_p$ ($p (x_p) =$ max($p$)) (dashed red line) and (ii) the wiping point $x^*$ (continuous white line). The streamwise coordinate is normalized with the jet opening $d$ and time axis is scaled as $\hat {t}= t / [t]$. The region delimited by these two points is referred to as wiping region, i.e. the region where most of the wiping work occurs. Here, the characteristic lines are gently curved because of the wave acceleration due to the strong gas–liquid interaction. These characteristic lines are continuous in case 1, while the waves occasionally reverse direction in case 3. This is probably due to the reduced effect of the jet flow outside the wiping region, combined with a temporary win of the viscous drag over gravity and inertia. In all cases, it appears that the amplitude of the jet oscillation (linked to the jet structure ‘$d$’) has a minor role compared with the large fluctuations experienced by the wiping point (linked to the jet structure ‘$v_r$’). These fluctuations are characterized by a ‘slow’ downward dynamic and a ‘fast’ upward dynamic, as if the mechanisms governing these two stages were completely different. The downward shift of the wiping point occurs at the scale of the wave advection time in the falling liquid waves, while the upward shift occurs at the scale of the gas jet advection time. It is noteworthy that the waves are propagated at approximately the substrate speed (${\hat {u} \approx 1}$) in both the final film and run-back flow, regardless of the wiping strength ($\varPi _g$), film intrusiveness ($h_0 / Z$) and characteristic speed of the side jet ($U_{sj}$) – which can reach values up to $\hat {u}_{sj} = 20$ in our simulations. This observation is in agreement with Lunz & Howell (2018), who investigated a falling liquid film subjected to external pressure perturbations, and showed that the maximum amplification occurs when the pressure disturbance moves at a speed comparable to the wave propagation.

Figure 13. Spatiotemporal contours of the normalized film thickness $\check {h}$ for cases (a) 1 and (b) 3, complemented with the temporal evolution of the impingement point $x_p$ (dashed red line) and the wiping point $x^*$ (continuous white line), both defined in the text. (a) Case 1 (DG); (b) case 3 (W).

From this perspective, the strip velocity plays a dominant role, while the wiping strength and the jet velocity plays a secondary role in selecting the coupling frequency. Tracking the white lines in figure 13, it appears that the time scale of the coupling is associated with (i) the time required to produce a wave in the run back flow, (ii) the time required for this wave to exit the wiping region. Below this zone, the wiping point is no longer able to ‘stick’ onto the wave's back and so it recovers its upward position, where a new wave will be generated. This mechanism is consistent with experimental data and the qualitative description in Mendez et al. (2019b), and justifies why a purely liquid-based scaling in which $[u]=U_p$ yields similar dimensionless frequencies in largely different wiping regimes.

Finally, to conclude on the similarity between the investigated conditions, we analyse the time response of the liquid film thickness to the unsteadiness of the wiping point for both cases 1 and 3. Figure 14 compares the evolution of the emPOD filtered film thickness and pressure gradient at the time averaged wiping point $\bar {x}^*$. The latter is defined as the location at which the time-averaged pressure gradient is maximum, i.e. $\bar {x}^*$ ($\overline {\partial _x p})(\bar {x}^*)=\mbox {max}(\overline {\partial _x p})$. The similarity between the two plots suggests that a comparable coupling mechanism takes place in the two cases despite the large differences in $\textit {Ka}$, wiping conditions and intrusiveness of the film. Because the dimensionless frequencies of the waves remain within the same range, we may infer that the interaction is mostly driven by the liquid film. Considering that the dimensionless groups related to the liquid film dynamics are in fair similarity between water and zinc (table 2), these results suggest that the coupling dynamics observed in water might not differ significantly from the one occurring in galvanization. It is striking that the wavy defects observed on galvanized products have typical wavelengths in the range 10–15 mm, which corresponds to dimensionless frequencies of $\hat {f} = (0.11\unicode{x2013}0.165)$, totally in line with what is found here for largely different conditions.

Figure 14. Lissajous curves for the dimensionless pressure gradient $\partial _{\hat {x}} \hat {p}$ and thickness $\hat {h}$ at the time averaged wiping point $\bar {x}_*$ for case (a) 1 and (b) 3. (a) Case 1 (DG); (b) case 3 (W).

5. Conclusions

We have numerically investigated the two-phase coupling instability taking place between an impinging gas jet and a liquid film dragged by an upward moving substrate. The two-phase simulations with a high Kapitza liquid allowed for characterizing the gas–liquid interaction in an unexplored wiping regime and much closer to the industrial conditions, especially from the liquid perspective.

The results show a strong gas–liquid interaction which produces undulations on the final film, both with low and high Kapitza liquids. The waves emerge 2-D in the impingement region, even when the run-back flow is 3-D. It suggests that the 2-D–3-D transition of the liquid waves is due to an intrinsic instability typical of a falling film at these Reynolds numbers, and not due to the gas jet, which is fundamentally 2-D. The correlation analysis in the gas jet reveals two flow structures acting at the time scale of the leading wave pattern: a symmetric oscillation around the jet axis (pattern ‘$d$’), and a deflection of the lower side jet triggered by the periodic formation of waves (pattern ‘$v_r$’). It is shown that, in all cases, the second mechanism has a stronger impact on the stresses at the interface and, thus, on the formation of the waves.

Remarkably, it is found that the velocity of the waves is similar to the substrate speed in both the final film and the run-back flow regardless of the wiping conditions, i.e. film intrusiveness, wiping strength and characteristic speed of the side jet. This observation suggests that the liquid film plays a dominant role in the coupling, and justifies the ‘appropriateness’ of a purely liquid based scaling that uses the substrate speed as a reference. On the other hand, the gas jet adjusts the characteristic frequency of the mechanism through the wiping strength, which modifies the time required to develop a run-back wave. From this perspective, the gas jet structures ‘$v_r$’ and ‘$d$’ can be seen as a consequence rather than a cause in the coupling.

In spite of the very different flow regimes analysed in this work, it is remarkable that the dynamics of the gas–liquid interaction is qualitatively similar. Although the system locks at a certain frequency that depends on both the gas jet and liquid film, the present results provide new insights about the scaling laws and the physics of the gas–liquid coupling.

Additionally, the data from the simulations allowed for assessing the validity of liquid film models in the context of high film Reynolds numbers and wiping strengths. The steady-state solutions of the liquid-film models compare well with the time-averaged profile from the CFD, and the assumptions underpinning the models are well justified. These findings open a promising path towards the application of liquid-film models in jet wiping, since it is currently the only possible alternative to simulate industrial galvanization with zinc.