1. Introduction
Integral boundary layer models for falling liquid films have been extensively used to study flow configurations that are encountered in many coating, chemical, heat and mass transfer processes. These models, also referred to as low-dimensional models, reduce the number of variables governing the problem by eliminating the velocity and pressure fields from the full set of Navier–Stokes equations, thus describing the dynamics of the liquid film as a function of thickness and flow rate. From the pioneering two-dimensional formulations proposed by Kapitza (Reference Kapitza1948a,b Reference Kapitza), Shkadov (Reference Shkadov1971) and extended by Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000), to the three-dimensional models firstly proposed by Demekhin & Shkadov (Reference Demekhin and Shkadov1985) and extended by Scheid, Ruyer-Quil & Manneville (Reference Scheid, Ruyer-Quil and Manneville2006), the literature on the topic is vast and discussed in various monographs (Alekseenko, Nakoryakov & Pokusaev Reference Alekseenko, Nakoryakov and Pokusaev1994; Chang & Demekhin Reference Chang and Demekhin2002; Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) and reviews (Chang & Demekhin Reference Chang and Demekhin1996; Craster & Matar Reference Craster and Matar2009; Ruyer-Quil et al. Reference Ruyer-Quil, Kofman, Chasseur and Mergui2014).
The capability of low-dimensional models to describe the dynamics of the liquid interface has been largely demonstrated for the fundamental case of a gravity-driven isothermal film (Alekseenko, Nakoryakov & Pokusaev Reference Alekseenko, Nakoryakov and Pokusaev1985; Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000, Reference Ruyer-Quil and Manneville2002; Scheid et al. Reference Scheid, Ruyer-Quil and Manneville2006), for which many experimental (Liu & Gollub Reference Liu and Gollub1994; Liu, Schneider & Gollub Reference Liu, Schneider and Gollub1995; Alekseenko, Markovich & Shtork Reference Alekseenko, Markovich and Shtork1996; Nosoko et al. Reference Nosoko, Yoshimura, Nagata and Oyakawa1996; Tihon et al. Reference Tihon, Serifi, Argyriadi and Bontozoglou2006; Dietze, Al-Sibai & Kneer Reference Dietze, Al-Sibai and Kneer2009; Mendez, Scheid & Buchlin Reference Mendez, Scheid and Buchlin2017b) and numerical (Salamon, Armstrong & Brown Reference Salamon, Armstrong and Brown1994; Gao, Morley & Dhir Reference Gao, Morley and Dhir2003; Nosoko & Miyara Reference Nosoko and Miyara2004; Dietze, Leefken & Kneer Reference Dietze, Leefken and Kneer2008; Malamataris & Balakotaiah Reference Malamataris and Balakotaiah2008; Meza & Balakotaiah Reference Meza and Balakotaiah2008; Doro & Aidun Reference Doro and Aidun2013; Dietze et al. Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014) investigations have been carried out. Because of their minor computational cost, if compared with full simulations, and because of the analytic insights they enable, these models have been largely used in more complex configurations. These include, among others, liquid films in the presence of interface shear stress (Frank Reference Frank2006, Reference Frank2008; Gatapova & Kabov Reference Gatapova and Kabov2008; Vellingiri et al. Reference Vellingiri, Tseluiko, Savva and Kalliadasis2013; Samanta Reference Samanta2014; Lavalle et al. Reference Lavalle, Vila, Lucquiaud and Valluri2017), for which the linear stability analysis based on the full Orr–Sommerfeld equations is presented by Lavalle et al. (Reference Lavalle, Li, Mergui, Grenier and Dietze2019).
Moreover, despite the restrictive hypotheses in their derivation, low-order formulations have been successfully validated with experimental and numerical data in operating conditions that are well outside their theoretical range of validity (Denner et al. Reference Denner, Charogiannis, Pradas, Markides, van Wachem and Kalliadasis2018). This has made integral model reliable tools to explore complex phenomena such as the origin of capillary ripples (Dietze Reference Dietze2016), the onset of circulating waves and flow reversal in inclined liquid films (Rohlfs & Scheid Reference Rohlfs and Scheid2014; Rohlfs, Pischke & Scheid Reference Rohlfs, Pischke and Scheid2017), the effect of co-flowing turbulent gas on the interface dynamics (Vellingiri et al. Reference Vellingiri, Tseluiko, Savva and Kalliadasis2013) or the formulation of active feedback flow control methods to suppress interface instabilities (Thompson, Tseluiko & Papageorgiou Reference Thompson, Tseluiko and Papageorgiou2015; Thompson et al. Reference Thompson, Gomes, Pavliotis and Papageorgiou2016; Tomlin et al. Reference Tomlin, Gomes, Pavliotis and Papageorgiou2019).
This work extends the classical integral boundary layer models for liquid films to a flow configuration for which these have never been used: the jet wiping process. This process consists in using an impinging gas jet to control the thickness of a coating film on a vertical moving substrate, and it is characterized by an unstable dynamics (Mendez et al. Reference Mendez, Gosset, Myrillas and Buchlin2017a), recently investigated experimentally by Gosset, Mendez & Buchlin (Reference Gosset, Mendez and Buchlin2019) and Mendez, Gosset & Buchlin (Reference Mendez, Gosset and Buchlin2019). In particular, it has been shown that instabilities on the gas jet can propagate to the impinged liquid and produce a non-uniform coating distribution referred to as undulation. Although several working hypotheses have been proposed, the mechanisms through which unsteadiness in the jet propagates to the liquid film are still not fully understood and are explored in this work.
The modelling of this configuration presents two distinctive features that have not been considered in the liquid film modelling literature: (1) the upward motion of the vertical substrate and (2) the simultaneous presence of time-dependent sources of shear stress and pressure gradient. Simplified theoretical modelling of the jet wiping have been proposed by Thornton & Graff (Reference Thornton and Graff1976), Tuck (Reference Tuck1983), Tuck & Vanden-Broeck (Reference Tuck and Vanden-Broeck1984), Tu & Ellen (Reference Tu and Ellen1986), Ellen & Tu (Reference Ellen and Tu1984), and Buchlin (Reference Buchlin1997); experimental and numerical validations have been provided by Lacanette et al. (Reference Lacanette, Gosset, Vincent, Buchlin and Arquis2006) and Gosset & Buchlin (Reference Gosset and Buchlin2007). These first formulations aimed at describing the mean thickness distribution of the liquid film under the action of the pressure gradient and the shear stress produced by the jet impingement, and, thus, at predicting the final coating thickness as a function of all the operating parameters.
The first works discussing the stability of the problem have been presented by Ellen & Tu (Reference Ellen and Tu1983) and Tuck (Reference Tuck1983). The first presented a linear stability analysis; the second discussed the possible evolution of kinematic waves on the liquid coat. Since then, most of the investigations on the process have been based on high fidelity numerical simulations, combining large eddy simulation (LES) of the gas jet with volume of fluid (VOF) treatment of the liquid film (Myrillas et al. Reference Myrillas, Gosset, Rambaud and Buchlin2009, Reference Myrillas, Rambaud, Mataigne, Gardin, Vincent and Buchlin2013; Eßl et al. Reference Eßl, Pfeiler, Reiss, Ecker, Riener and Angeli2017; Pfeiler et al. Reference Pfeiler, Eßl, Reiss, Riener, Angeli and Kharicha2017; Aniszewski et al. Reference Aniszewski, Zaleski, Popinet and Saade2019). While these simulations can potentially provide a complete picture of the unstable interaction between the jet and the gas flow, their computational cost remains prohibitively large for analysing configurations of industrial interest, as discussed by Aniszewski et al. (Reference Aniszewski, Zaleski, Popinet and Saade2019).
A theoretical analysis of the stability of the process has been proposed by Hocking et al. (Reference Hocking, Sweatman, Fitt and Breward2010), who used a quasi-steady formulation to study the evolution of liquid film disturbances. Hocking and coworkers concluded that the coating film is neutrally stable and incapable of producing the undulation patterns observed in the wiping lines without the presence of disturbances produced by the gas jet. Using the same quasi-steady formulation, Johnstone et al. (Reference Johnstone, Kosasih, Phan, Dixon and Renshaw2019) investigated the response of the liquid coat to a set of possible unsteady behaviour of the impinging jet. The formulation presented in these works neglects the role of inertia in the liquid film, disregarding the nonlinear contribution of advection.
In this work the extension of more advanced integral models to the jet wiping process is used to study the dynamic response of the liquid film to various disturbances on the gas jet. These include localized perturbation, simulated by pulsation of the wiping actuators, and various kinds of harmonic and non-harmonic oscillations. The general form of these models is presented in § 2, while §§ 3 and 4 provide the details of the laminar and turbulent models. Appendix A provides complementary material for the full derivation of the models.
Among the laminar models in § 3, this work covers the zero-order (ZO) formulation of the jet wiping, the extension of the integral boundary layer (IBL) model from Kapitza–Shkadov (Shkadov Reference Shkadov1971; Shkadov & Beloglazkin Reference Shkadov and Beloglazkin2017) and the extension of the weighted integral boundary layer (WIBL) from Ruyer-Quil and Manneville (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000, Reference Ruyer-Quil and Manneville2002). The proposed transition and turbulence model in § 4 combines ideas from mixing length formulation (van Driest Reference van Driest1956; King Reference King1966; Geshev Reference Geshev2014) and shallow water formulations (James et al. Reference James, Lagrée, Le and Legrand2019; De Vita et al. Reference De Vita, Lagrée, Chibbaro and Popinet2020).
Section 5 reviews the implementation of the wiping actuators in the models. Section 6 presents the numerical methods, including the finite volume (FV) solver implemented to validate the integral models (in § 6.1) and the direct numerical simulations using OpenFoam (in § 6.2) that were implemented to set-up simple validation test cases. The results are presented in § 7, including the numerical validation (§ 7.1), the relative weight of all the forces governing the process (§ 7.2) and the frequency response of the coating film (§ 7.3). Finally, the impact of the modelling strategy on the identified wave formation mechanisms is discussed in § 7.4. Conclusions and perspectives are presented in § 8.
2. The integral formulation for the jet wiping
The jet wiping process is represented schematically in figure 1. A liquid film is dragged along a vertical plate moving upward at a constant velocity 
$U_p$ and is impinged upon by a gas jet. The configuration is assumed two dimensional, with incompressible liquid flow bounded by the plate at 
$y=0$, and the dynamic liquid interface at 
$y=h(x,t)$.
Figure 1. Schematic of the jet wiping process: a nozzle with an opening 
$d$ and a stagnation pressure 
${\rm \Delta} P_N$ releases a jet flow at a distance 
$Z$ from a dip-coated substrate moving at a speed 
$U_p$. The impingement produces a wiping meniscus (region WR). This forces a run-back film to flow backward (region RB) and leaves a thinner liquid film downstream (region FC) before solidification takes place.
The origin 
$x=0$ is located at the nozzle axis, and the streamwise coordinate 
$x$ is oriented in the direction of gravity and counter to the substrate velocity. The impinging jet flow produces a pressure 
$p_g(x,t)$ and a shear stress distribution 
$\tau _g(x,t)$ that identify three areas, qualitatively pictured in figure 1 on the right. In the wiping region (indicated as WR), the pressure gradient imposed by the gas jet forces part of the liquid to reverse direction resulting in a wiping meniscus. The falling liquid forms the run-back flow in the region 
$x\rightarrow \infty$ (indicated as RB); the remaining liquid evolves upward in the final coating region 
$x\rightarrow -\infty$ (indicated as FC).
In a one-way coupling formulation, it is assumed that the presence of the liquid film does not influence the gas jet. This assumption has been extensively validated for the prediction of the averaged final coating thickness (Lacanette et al. Reference Lacanette, Gosset, Vincent, Buchlin and Arquis2006; Gosset & Buchlin Reference Gosset and Buchlin2007), but it is certainly not able to simulate the complex interaction between the two flows analysed by Gosset et al. (Reference Gosset, Mendez and Buchlin2019) and Mendez et al. (Reference Mendez, Gosset and Buchlin2019). In this work, this formulation is used to de-couple the dynamics of the liquid from the one of the gas jet and to analyse the liquid film frequency response and possible mechanisms of undulation formation. We thus assume that both the pressure and the shear stress produced by the jet solely depend on the nozzle gauge stagnation pressure 
${\rm \Delta} P_N$, the nozzle opening 
$d$, its discharge coefficient 
$C_d$ and the standoff distance 
$Z$. This dependency is described in § 5.
All the integral models investigated in this work rely on the Navier–Stokes equation and the related boundary conditions in the ‘long-wave’ formulation. This formulation is derived by scaling the cross-streamwise direction with a reference length 
$[h]$, which is much smaller than the streamwise reference length 
$[x]$. Appendix A summarizes this derivation and presents the rationale behind the choice of the reference quantities used to scale the problem. For conciseness, these reference quantities are listed in table 1, and the focus is here kept on the derivation of the various integral models. In what follows, dimensionless variables scaled with respect to the quantities in table 1 are indicated with a hat (e.g. 
$\hat {h}=h/[h]$).
Table 1. Reference quantities for the Shkadov-like scaling, for which 
$\varepsilon =Ca^{1/3}$, with 
$Ca=\mu _l U_p/\sigma$ the capillary number.
In the long-wave formulation of the problem, the dimensionless continuity and momentum equations in the 
$x$ and the 
$y$ directions reduce to the boundary layer equations:
\begin{gather} \partial_{\hat{x}} \hat{u}+ \partial_{\hat y} \hat{v}=0, \end{gather}
\begin{gather}\varepsilon Re(\partial_{\hat t}\hat{u}+\hat {u} \partial_{\hat x} \hat{u} + \hat {v} \partial_{\hat y} \hat{u}) =- \partial_{\hat x} \hat{p}_l+\partial_{\hat y \hat y} \hat{u} +1, \end{gather}
\begin{gather}0= \partial_{\hat y } \hat{p}_l. \end{gather}
Here 
$\hat {p}_l$ is the pressure in the liquid, 
$\hat {u}$ and 
$\hat {v}$ are the streamwise and cross-stream velocity components, 
$\varepsilon =[h]/[x]=Ca^{1/3}$ is the film parameter, with 
$Ca=\mu _l U_p/\sigma$ the capillary number and 
$Re=[u][h]/\nu _l=(U_p^3/g\nu _l)^{1/2}$ the global Reynolds number of the process, to be distinguished from other Reynolds numbers that will be introduced later.
Since the proposed scaling laws hold for 
$\varepsilon \ll 1$, the long-wave formulation presented in this work is valid for 
$Ca^{1/3}\ll 1$. This generally occurs in galvanizing conditions, where typically 
$\mu _l\approx 0.003\,\textrm{Pa s}$, 
$\sigma \approx 0.8$ N m
$^{-1}$ and 
$U_p=1{-}2$ m s
$^{-1}$ and, hence, 
$Ca\approx 0.004{-}0.008$. The formulation proposed in this work is also valid for the experimental conditions encountered in the Essor laboratory (see Buchlin Reference Buchlin1997) developed at the von Karman Institute (VKI), operating with water (
$\mu _l\approx 0.001$ 
$\textrm{Pa}\ \textrm {s}$, 
$\sigma \approx 0.07$ N m
$^{-1}$ at 
$U_p=0.2{-}2$ m s
$^{-1}$, i.e. 
$Ca=0.003{-}0.03$). On the other hand, the experiment recently conducted in the VKI Ondule laboratory (Gosset et al. Reference Gosset, Mendez and Buchlin2019; Mendez et al. Reference Mendez, Gosset and Buchlin2019) using dipropilene glycole (
$\mu _l \approx 0.1\,\textrm {Pa} \ \textrm {s}$, 
$\sigma \approx 0.03$ N m
$^{-1}$ at 
$U_p=0.2{-}0.4$ m s
$^{-1}$, i.e. 
$Ca\approx 0.6{-}1.3$) requires a different scaling strategy. While the previous experimental work of the authors has focused on the dynamics of very viscous flows (see also Mendez et al. Reference Mendez, Scheid and Buchlin2017b), this work focuses on the low 
$Ca$ limit, in which surface tension plays a more important role.
The kinematic boundary conditions at the wall and at the gas–liquid interface set:
\begin{equation} \begin{cases} \hat{\boldsymbol{v}}=(\hat u,\hat v)=(-1,0) & \mbox{in } \hat{y}=0, \\ \hat v= \partial_{\hat t} {h}+\hat{u} \partial_{\hat x} \hat h & \mbox{in } \hat{y}=\hat{h}. \end{cases} \end{equation}The dynamic boundary conditions (see (A3) and (A9)) at the interface simplify to
\begin{gather} \hat p_l|_{\hat h}=\hat p_g(\hat x,\hat t)-\partial_{\hat x \hat x}\hat h \mbox{ in } \hat{y}=\hat{h}, \end{gather}
\begin{gather}\partial_{y} \hat u|_{\hat h}=\hat{\tau}_g (\hat x,\hat t) \mbox{ in } \hat{y}=\hat{h}. \end{gather}All integral models reduce the modelling complexity by rendering the problem one dimensional. Integrating (2.1a)–(2.1c) across the film thickness using Leibniz integral rule together with the boundary conditions (2.3a) and (2.3b) gives
\begin{gather} \partial_{\hat t} \hat h+\partial_{\hat x }\hat q=0, \end{gather}
\begin{gather}\varepsilon Re (\partial_{\hat t} \hat q+ \partial_{\hat x} {\mathcal{F}})=\hat h(1-\partial_{\hat x} \hat{p}_g+\partial_{\hat x\hat x\hat x}\hat h) +{{\rm \Delta} \hat \tau}, \end{gather}
where 
$\hat q$ is the volumetric flow rate per unit width, and
\begin{equation} \mathcal{F}=\int_0^{{\hat h}}\hat{u}^2\,d{\hat y} \quad \mbox{and} \quad{\rm \Delta} \hat \tau \equiv \hat{\tau}_g-\hat{\tau}_w=\hat{\tau}_g-\partial_{\hat y} \hat u|_{\hat y=0} \end{equation}
are the advection and the shear stress terms, respectively, with 
$\hat {\tau }_w$ the wall shear stress. To determine the functional forms of these (and, thus, to close integral models), some assumptions on the velocity profile are required.
In this work we assume that the velocity profile is the superposition of three terms:
\begin{equation} \hat{u}(\hat{x},\hat{y},\hat{t})=\hat{u}_{F}(\hat{x},\hat{y}, \hat{t})+\hat{u}_{C}(\hat{x},\hat{y},\hat{t})+\hat{u}_{P}= \hat{u}_{F}(\hat{x},\hat{y},\hat{t})+\hat{\tau}_g(\hat{x},\hat{t})\hat y -1. \end{equation}
The term 
$\hat {u}_{F}$ accounts for the contributions of gravity, viscous stresses, surface tension and the pressure gradient. The term 
$\hat {u}_{C}=\hat {\tau }_g\hat y$ accounts for the shear stress produced at the gas liquid interface, while 
$\hat {u}_{P}=-1$ accounts for the motion of the substrate. This kinematic decomposition satisfies the boundary conditions if 
$\hat {u}_{F}=0$ at 
$\hat {y}=0$ and 
$\partial _{\hat {y}}\hat {u}_{F}=0$ at 
$\hat {y}=\hat {h}$. For later convenience, it is interesting to identify the reference velocity and the associated Reynolds number for each of the three contributions. The flow rate per unit width can be split accordingly as
\begin{equation} \hat q\equiv\int_0^{\hat h} \hat u \, \textrm{d} \hat y=\hat{q}_F+\hat{q}_C+\hat{q}_P=\hat{q}_F+\frac{1}{2}{\hat{\tau}_g \hat {h}^2}-\hat{h}, \end{equation}
from which the associated local (i.e. function of 
$\hat {x}$) Reynolds numbers are
\begin{equation} Re_F=\frac{q_F}{\nu}=|\hat{q}_F| Re;\quad Re_{\tau}=\frac{1}{2}{\hat h^2 |\hat{\tau}_g}| Re;\quad Re_h=\hat{h} Re . \end{equation} It is worth noting that the term 
$\hat {q}_F$ corresponds to a falling film flow in the absence of the other terms, but it can eventually lead to a negative contribution 
$\hat {q}_F<0$ in a strongly dominated shear stress flow (if 
$\hat {h}^2 |\hat {\tau }_g|\gg |\hat {q}_F|$), as is the case in the run-back flow region (more about this in § 4.1). The models developed in the following sections only differ in the treatment of this term.
3. Laminar film models
The term 
$\hat {u}_{F}$ in (2.6) is decomposed in a series of basis functions as
\begin{equation} \hat{u}_{F}=\sum^{N}_{j=0} a_j(\hat{x},\hat{t}) f_j\left(\frac{\hat y}{\hat h(\hat x,\hat t)}\right), \end{equation}
where 
$a_0$ is of 
${O}(1)$ and 
$a_j$ with 
$j>0$ are corrections of 
${O}(\varepsilon )$. Following Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2000, Reference Ruyer-Quil and Manneville2002) and Ruyer-Quil et al. (Reference Ruyer-Quil, Kofman, Chasseur and Mergui2014), the basis functions 
$f_j$ are taken as
\begin{equation} f_j=\bar{y}^{\,j+1}-\frac{j+1}{j+2} \bar{y}^{\,j+2}, \end{equation}
where 
$\bar {y}=\hat {y}/\hat {h}(x,t)$ is the reduced coordinate. This choice was introduced for falling liquid films, imposing that each basis function satisfies the boundary conditions. This enables reduced-order models based on Galerkin projections (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). The flow rate per unit width, highlighting the contributions in (2.7), becomes
\begin{equation} \hat q_F\equiv\int_0^{\hat h} \hat u_F \, \textrm{d} \hat y=\sum^{N}_{j=0} \frac{2}{(j+2)(j+3)}a_j. \end{equation} In all the laminar models presented in this work, valid at 
${O}(\varepsilon )$, only the first term (
$\,j=0$) contributes to the advection term 
$\mathcal {F}$, since this is multiplied by 
$\varepsilon$ in (2.4). The advection term, using (2.6) and (3.1), reads as
\begin{equation} \hat{\mathcal{F}}\equiv\int^{\hat{h}}_{0} \hat{u}^{2}\, \textrm{d} \hat y =\frac{1}{3} \hat{\tau}^{2}_g \hat h^3+\frac{5}{12} \hat{\tau}_g a_0 \hat {h}^2 -{\hat{\tau}_g} h^2 +\frac{2}{15} h{a^{2}_0}-\frac{2 }{3}a_0 \hat h+\hat h. \end{equation}
The shear stress contribution is solely linked to the first coefficient 
$a_0$ regardless of the number of terms included in the expansion of the velocity profile. The shear stress term, using (2.6) and (3.1), reads as
\begin{equation} {\rm \Delta} \hat \tau\equiv\hat{\tau}_g-\partial_{\hat y} \hat u|_{\hat y=0}=-\hat{\tau}_{wF}=-\frac{a_0 }{\hat h }, \end{equation}
where 
$\hat {\tau }_{wF}$ is the wall shear stress produced by the 
$\hat {u}_F$ portion of the velocity profile.
Before presenting the derivation of the complete WIBL model at 
${O}(\varepsilon )$, it is worth introducing the ZO formulation and the IBL formulation.
3.1. Zero-order (inertialess) formulation
This model is based on two assumptions. First, only the first 
$j=0$ term of the velocity profile is relevant, that is, 
$N=0$ in (3.1). Second, the inertial effects can be neglected, that is, 
$\varepsilon Re\sim 0$: the left-hand side of (2.4b) vanishes. The velocity profile is parabolic, i.e.
\begin{equation} \hat u=a_0(\hat x,\hat t) \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}
and the coefficient 
$a_0$ can be derived from the flow rate definition
\begin{equation} \hat q\equiv\int_0^{ \hat h} \hat u\,\textrm{d} \hat y=\frac{a_0 \hat h}{3}+\frac{ \hat h^2 \hat \tau_g}{2}- \hat h \rightarrow a_0=\frac{3 \hat{q}_F}{\hat h}= \frac{3 \hat q}{\hat h}-\frac{3}{2} \hat h \hat \tau_g+3. \end{equation}Using (3.7) in (3.5), the shear stress term becomes
\begin{equation} {\rm \Delta} \hat \tau\equiv -\hat{\tau}_{wF}=\frac{3}{2} \hat\tau_g-\frac{3 \hat q}{ \hat h^2}-\frac{3}{ \hat h}. \end{equation}Introducing this into (2.4b), and recalling that the left-hand side is set to zero, the flow rate is
\begin{equation} \hat q=\frac{ \hat h^3}{3} (1-\partial_{ \hat x}{\hat p}_g+\partial_{ \hat x \hat x \hat x} \hat h )+\frac{1}{2} \hat \tau_g \hat h^2 -\hat h. \end{equation}Introducing (3.9) in (2.4a) yields a single equation governing the film dynamics:
\begin{equation} \partial_{ \hat t} \hat h+\partial_{\hat x}\left[\frac{\hat h^3}{3} (1-\partial_{\hat x} \hat p_g+\partial_{\hat x\hat x\hat x} \hat h )+\frac{1}{2}\hat \tau_g \hat h^2 -\hat h\right]=0. \end{equation}This ZO model has been widely used in the literature of the jet wiping for linear stability analysis (Tuck Reference Tuck1983; Tu & Ellen Reference Tu and Ellen1986; Gosset Reference Gosset2007) or sensitivity studies similar to those performed in this work: Hocking et al. (Reference Hocking, Sweatman, Fitt and Breward2010) used this formulation to study the evolution of liquid disturbances for an ideally stationary jet; Johnstone et al. (Reference Johnstone, Kosasih, Phan, Dixon and Renshaw2019) used it to study the response of the liquid film to an oscillating jet.
Since this work focuses on more advanced formulations, the time-dependent simulation of the ZO model is not investigated further. It is nevertheless interesting to use this model to illustrate the basic features of a steady-state solution and the propagation of small flow disturbances, for which one could expect inertia to play a negligible role. In steady-state conditions (
$\partial _{ \hat t}\hat h=-\partial _{\hat x}\hat q=0$), and by neglecting the contribution of the surface tension term 
$\partial _{ \hat x \hat x \hat x} \hat h=0$ (which is known to have little influence on the final thickness in case of strong wiping, as shown by Tuck & Vanden-Broeck (Reference Tuck and Vanden-Broeck1984), Yoneda (Reference Yoneda1993) and Buchlin (Reference Buchlin1997)), (3.9) reduces to a cubic polynomial in 
$\hat {h}$.
At each location 
$\hat x$ (hence given 
$\partial _{\hat {x}} \hat {p}_g (\hat {x}),\hat {\tau }_g(\hat {x})$), and for a given flow rate 
$\hat {q}<0$, this polynomial admits a negative solution (of no interest) and two positive solutions 
$\hat h^{+}( \hat x)$ and 
$\hat h^{-}( \hat x)$. These branches of the positive solution give the liquid thickness in the final coating region 
$\hat h^{-}( \hat x)\rightarrow h_f$ for 
$\hat x \rightarrow -\infty$ and the run-back region 
$\hat h^{+}( \hat x)\rightarrow h_R$ for 
$\hat x \rightarrow \infty$. The admissible flow rate 
$\hat {q}$ can thus be computed by imposing that the two branches of solutions meet (see Tuck & Vanden-Broeck Reference Tuck and Vanden-Broeck1984; Hocking et al. Reference Hocking, Sweatman, Fitt and Breward2010) at a critical point 
$\hat {x}=\hat {x}_c$ to form a continuous thickness profile:
\begin{equation} \hat h(\hat{x})= \begin{cases} \hat{h}^{+}(\hat{x}) & \mbox{for } \hat{x}\leq\hat{x}_c, \\ \hat{h}^{-}(\hat{x}) & \mbox{for } \hat{x}>\hat{x}_c. \end{cases} \end{equation} A simple estimation of the final thickness can be obtained under the assumption that the system operates in optimal conditions (that is, 
$\partial _{ \hat h}{\hat q}=0$), assuming that the maximum pressure gradient and shear stress act at the same location 
$\hat x^{*}$. From (3.9), this gives
\begin{equation} \partial_{ \hat h}{\hat q}(\hat {h}^{*})=\hat {h}^{*2}(1-\partial_{\hat x}\hat p_{g}^{*})+\hat \tau_g^{*} \hat h^{*}-1=0, \end{equation}and, thus, the film thickness in this location is
\begin{equation} \hat h^{*}=\frac{-\hat{\tau}_g^{*}+\sqrt{{\hat \tau_g^{*2}}+4(1-\partial_{ \hat x}\hat{p}_{g}^{*})}}{2(1-\partial_{\hat x}\hat{p}_{g}^{*})}. \end{equation} It is worth noting that the downward orientation of the 
$\hat x$ axis in this work is opposite to the one used in the jet wiping literature (Tuck & Vanden-Broeck Reference Tuck and Vanden-Broeck1984; Gosset & Buchlin Reference Gosset and Buchlin2007; Hocking et al. Reference Hocking, Sweatman, Fitt and Breward2010), but in line with the falling liquid film literature (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012; Ruyer-Quil et al. Reference Ruyer-Quil, Kofman, Chasseur and Mergui2014). The wiping thus occurs in a region in which 
$\partial _{ \hat x}\hat p_{g}^{*}<0$ and 
$\hat \tau _g^{*}>0$. Introducing this value of the film thickness in (3.9), taking 
$\partial _{\hat x}\hat p_{g}=\partial _{\hat x}\hat p_{g}^{*}$ and 
$\hat \tau _g=\hat \tau _g^{*}$, allows for estimating the withdrawn flow rate 
$q(\hat h^{*})$. The final coating thickness 
$\hat {h}_f$ can thus be estimated using again (3.9) in the far-field conditions (where 
$\hat \tau _g= \partial _{ \hat x}\hat p_{g}=0$). This approach, known as the zero-dimensional (0-D) knife model, was proposed by Buchlin (Reference Buchlin1997) and validated on several numerical and experimental works (Lacanette et al. Reference Lacanette, Gosset, Vincent, Buchlin and Arquis2006; Gosset & Buchlin Reference Gosset and Buchlin2007).
The polynomial in the far-field condition is shown in figure 2 (see also Snoeijer et al. Reference Snoeijer, Ziegler, Andreotti, Fermigier and Eggers2008). For a wiping condition yielding 
$\hat {q}({h}^*)=-0.1$, the corresponding final thickness (
$\hat {h}_f$) and run-back flow thickness (
$\hat {h}_R$) are shown. At the lowest limit of the flow rate (
$\hat {q}=-2/3$), only one solution is admissible for the film thickness (corresponding to 
$\hat {h}=1$ in the chosen scaling). This corresponds to the well-known limit for the drag-out problem (Deryagin & Levi Reference Deryagin and Levi1964; Rio & Boulogne Reference Rio and Boulogne2017) in the gravity dominated regime.
Figure 2. Flow rate versus thickness relations for a film flowing along a vertical moving wall, assuming that 
$\partial _{\hat x} \hat {p}_g=\hat {\tau }_g=0$ in (3.9) and neglecting the surface 
$\partial _{\hat x\hat x\hat x}\hat {h}$. For a given flow rate 
$\hat q(h^{*})<0$, the thin and thick solutions are approximations of 
$\hat h_f=\lim _{\hat x \to +\infty } \hat {h}(\hat {x})$ and 
$\hat h_R=\lim _{\hat x \to -\infty } \hat {h}(\hat {x})$, respectively. Using 
$h^{*}$ from (3.13) yields the 0-D knife model of the wiping process.
Equation (3.10) allows for important considerations on the propagation of small disturbances on a flat film over an upward-moving substrate. Neglecting the surface tension contribution, this equation simplifies to
\begin{equation} \partial_{\hat t} \hat h+( \hat h^2-1) \partial_{ \hat x} \hat h=0, \end{equation}
that is, a standard kinematic wave equation (Whitham Reference Whitham1999): disturbances on a film 
$\hat h<1$ (that is, in the final coat) propagate at a negative velocity (that is, upward) while disturbances on a film 
$\hat h>1$ (that is, in the run-back flow) propagate at a positive velocity (that is, downward).
3.2. Integral boundary layer formulation
As for the ZO model, this formulation considers only the first term in the velocity profile (3.6), but it does not assume that the left-hand side of (2.4b) vanishes. Using (3.6) with the coefficient 
$a_0$ from (3.7), the advection term 
$\mathcal {F}$ in (2.5a,b) yields
\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}. \end{equation}The shear stress term remains the same as that in (3.8). The resulting model is an extension to the jet wiping process of the classical self-similar integral models proposed by Haar (Reference Haar1965), Shkadov (Reference Shkadov1971) and Shkadov & Beloglazkin (Reference Shkadov and Beloglazkin2017). To the best of the authors’ knowledge, this integral model has never been used in the analysis of the jet wiping process.
3.3. Weighted integral boundary layer model
The full expansion in (3.1) is now considered, taking up to 
$N=4$ terms. This number can be derived based on the order of magnitude analysis (see Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012, § 6.6). The method of weighted residuals to derive the coefficients 
$\mathcal {A}=\{a_0,a_1,a_2,a_3,a_4\}$ in (3.1) for a falling liquid film was proposed by Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2002). Denoting the momentum equation (2.1c) as an operator 
$\mathcal {M}(\hat {u})=0$, it is possible to construct the residuals 
$\mathcal {R}_j$ from the projections 
$\mathcal {R}_j=\langle w_j, \mathcal {M}(\sum _j a_j f_j(\hat {y}))\rangle$, with 
$\langle \cdot , \cdot \rangle$ denoting the continuous inner product over the space of possible solutions 
$u\in \mathcal {S}$ and 
$w_j$ a set of weight functions. Setting all residuals to 
$\mathcal {R}_j=0$ leads to a system of equations for the 
$a_j$ coefficients.
For the modelling of a film falling along a fixed wall and passive atmosphere, Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) compare various weighted residual methods, differing by choice of weight functions 
$w_j$. Their comparison shows that all methods converge, if four 
$w_j$ are taken – Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012) also show that some methods converge ‘faster’ than others: for instance, the Galerkin approach with 
$w_j=f_j$ converges with only one residual, while a collocation approach, with 
$w_j=\delta _j$, with 
$\delta _j$ a set of equally spaced Dirac functions, requires four residuals – to a model that can be derived following a polynomial matching procedure that involves neither weights nor residuals.
We here solely focus on this second approach, although this requires considerably more algebra, for two reasons. First, because this allows us to derive an explicit equation for the coefficients 
$a_j$, while this is not needed in the weighted residual approach. Second, because this approach does not depend on the choice of the 
$w_j$ and allows a discussion on convergence to be avoided. Nevertheless, we still refer to the resulting model as WIBL, in line with the literature on falling liquid films.
Introducing the expansion of the velocity profile (3.1) in the momentum equation (2.1b) yields a polynomial in the reduced coordinate 
$\bar {y}$:
\begin{equation} \mathcal{P}\left(x,t,\bar{y}=\frac{y}{h}\right)= \sum^{N}_{j=0}\,P_j(\mathcal{A},x,t) \bar{y}^j=0. \end{equation}
Because this polynomial must be identically null, all the functions 
$P_j (\mathcal {A},x,t )$ must be null. Moreover, we note that the viscous term on the right-hand side of (2.1b) decreases the degree of the polynomials by two while the left-hand side increases it by two. Therefore, only the 
$N=0$ should be introduced on the left-hand side, and up to 
$N=4$ terms should be introduced on the right-hand side, recalling that 
$a_0\sim {O}(1)$ and 
$a_i\sim {O}(\varepsilon ) \forall \,i \in [1,4]$. The resulting polynomial in (3.16) is thus of order four.
Setting all of its coefficients to zero leads to a system of five equations of the form
\begin{equation} [P_j(\mathcal{A},x,t)]^{4}_{0}=0 \rightarrow \varGamma \mathcal{A}=\mathcal{G}, \end{equation}where the system matrix reads as
\begin{equation} \varGamma=\frac{1}{h^2} \begin{pmatrix} 1 & -2 & 0 & 0 & 0 \\ 0 & 4 & -6 & 0 & 0\\ 0 & 0 & 9 & -12 & 0 \\ 0 & 0 & 0 & 16 & -20\\ 0 & 0 & 0 & 0 & 25 \end{pmatrix}, \end{equation}
and the vector 
$\mathcal {G}=\{G_0,G_1,G_2,G_3,G_4\}$ is shown in table 5 of appendix B. The solution of the system leads to the full expression of the coefficient, given in table 6 of appendix B. All the coefficients, except 
$a_0$, are then introduced in (3.3) to identify the link between 
$a_0$ – thus, the wall shear stress term in (3.5) – and the flow rate. The resulting expression, following the classical notation from asymptotic expansions (Howison Reference Howison2005), is of the form
\begin{equation} \hat q=\hat q^{(0)} (a_0, \hat h)+\varepsilon Re \, q^{(1)}(\hat{h},a_0,\partial_{\hat t} a_0, \partial_{\hat x} a_0, \partial_{\hat x} \hat h, \partial_{\hat t} \hat h ), \end{equation}
having considered only the functional dependency on the unknown variables. This expression can be used to compute 
$a_0$ following the same notation:
\begin{equation} a_0=a_0^{(0)} (\hat q, \hat h)+\varepsilon Re \,a_0^{(1)}(\hat{h},a_0,\partial_{\hat t} a_0, \partial_{\hat x} a_0, \partial_{\hat x} \hat h, \partial_{\hat t} \hat h ). \end{equation} Both expressions (3.19) and (3.20) are shown in table 7 of appendix B in their complete forms ((B11) and (B12), respectively). As expected, the first-order term recovers the coefficient from the leading-order model in (3.7), while the other represents 
${O}(\varepsilon )$ corrections. At this stage, the coefficient 
$a_0$ is still implicitly defined. However, if 
$Re\ll 1/\varepsilon$, the asymptotic expansion framework allows for substituting 
$a_0\approx a_0^{(0)}+{O}(\varepsilon )$ in 
$a_0^{(1)}$ and neglecting higher-order terms, as done in (3.19) and (3.20). The resulting shear stress is
\begin{align} {\rm \Delta} \tau&=\frac{3}{2}\tau_g-\frac{3q}{h^2}-\frac{3}{h}+\varepsilon Re \left(-\frac{19 {{h }^{3}} {{\tau}_g} \partial_{x} {{\tau}_g}}{3360}-\frac{17 h q \partial_{x} {{\tau}_g} }{560}- \frac{3 {{h }^{2}} \partial_{x} {{\tau}_g}}{560}-\frac{ {{h }^{2}} \partial_{t} {{\tau}_g} }{40}\right.\nonumber\\ &\quad -\frac{ h {{\tau}_g} \partial_{x} q}{56}-\frac{18 q \partial_{x} q }{35 h }-\frac{4 \partial_{x} q }{35}-\frac{ \partial_{t} q }{5}-\frac{ {{h }^{2}} {{{{\tau}_g}}^{2}} \partial_{x} h}{112}-\frac{ q {{\tau}_g} \partial_{x} h }{280}-\frac{3 h {{\tau}_g} \partial_{x} h }{140}\nonumber\\ &\quad \left.+\frac{12 {{q }^{2}} \partial_{x} h}{35 {{h }^{2}}}+\frac{6 q \partial_{x} h }{35 h }+\frac{ \partial_{x} h }{35}\vphantom{\frac{19 {{h }^{3}} {{\tau}_g} \partial_{x} {{\tau}_g}}{3360}} \right). \end{align} The WIBL model for the jet wiping problem is obtained by introducing (3.21) and (3.15) into (2.4a) and (2.4b). Observe that the 
${O}(1)$ terms in (3.21) are those in the IBL model. Among the fourteen 
${O}(\varepsilon )$ terms, it is worth noticing the one involving the partial time derivative of the flow rate (the eighth term). For computational purposes, it is convenient moving this term on the left-hand side of (2.4b), resulting in a coefficient 
$\beta =6/5$ on the term 
$\partial _t \hat {q}$ (see § 6).
4. The transition and turbulent boundary layer
A falling liquid film is laminar for Reynolds numbers below 
$\approx 100$, and it is turbulent above 
$\approx 400$ (Ishigai et al. Reference Ishigai, Nakanisi, Koizumi and Oyabu1972; Alekseenko et al. Reference Alekseenko, Nakoryakov and Pokusaev1994; Karimi & Kawaji Reference Karimi and Kawaji1999). The intermediate range is the transition region, and cannot be described solely in terms of Reynolds number. Laminar integral boundary layer models have been proved successful (see Denner et al. Reference Denner, Charogiannis, Pradas, Markides, van Wachem and Kalliadasis2018 for cases up to 
$Re\approx 80$) well above their theoretical range of validity (which sets 
$Re\sim {O}(1)$), while much higher Reynolds numbers, like those considered in this work, might need a different treatment.
A classic theoretical formulation for turbulent liquid films is based on mixing length theory (van Driest Reference van Driest1956; King Reference King1966; Geshev Reference Geshev2014) and a Reynolds-averaged formulation of the velocity field in which the effect of turbulence is modelled by an additional eddy viscosity. This formulation is based on the statistically stationary assumption and is of difficult extension to the integral formulation of interest to this work. A different approach is commonly encountered in the literature of shallow water flows (see James et al. Reference James, Lagrée, Le and Legrand2019), in which the most celebrated empiricism consists in introducing a correlation for the wall shear stress (see also Katopodes Reference Katopodes2018), and a shape factor for the velocity profile. As these allow for keeping the integral nature of the model formulation, a similar approach is pursued in this work.
The transition and turbulent boundary layer (TTBL) model for the jet wiping problem proposed in this work makes no pretension of completeness; on the contrary, it offers a first attempt to analyse the possible impact of turbulence on the response of the coating thickness. Following the same self-similarity argument supporting the IBL model, the proposed model extends the IBL model to a turbulent liquid film. A similar extension of the WIBL for a falling liquid film, using the mixing length formulation, is presented by Mukhopadhyay, Chhay & Ruyer-Quil (Reference Mukhopadhyay, Chhay and Ruyer-Quil2017). The proposed closures for the wall shear stress and the advection terms are described in §§4.1 and 4.2, respectively.
4.1. Closure for the wall shear stress
Following the shallow water literature, the wall friction is modelled in terms of a friction coefficient 
$C_f$ that is a function of the (local) Reynolds number. In the jet wiping problem, where the liquid streamwise velocity component is composed of the terms in (2.6), one must first establish which of the Reynolds numbers in (2.8a–c) controls the transition to turbulence and the flow regime. Moreover, since the thickness of the liquid film varies significantly between the final coating region (
$\hat {h}\ll 1$) and the run-back flow region (
$\hat {h}\sim 1$), different regimes (laminar, transition, fully turbulent) can be expected in different regions. We begin by introducing the skin friction coefficient in the laminar models. From (3.6) and (3.7), the wall shear stress reads as
\begin{equation} { \partial_{\hat{y}}\hat{u}|_{\hat {y}=0}=\frac{3 \hat{q}_F}{\hat{h}^2}+\hat{\tau}_g=\hat{\tau}_{wF}(\hat{q}_F)+\hat{\tau}_g}. \end{equation} Introducing the skin friction coefficient based on the mean velocity 
$\hat {q}_F/\hat {h}$, we obtain the dimensional and dimensionless friction terms 
$\hat {\tau }_{wF}$,
\begin{equation} { {\tau}_{wF}=\frac{1}{2}\rho \frac{q_F|q_F|}{h^2} C_f\Longleftrightarrow \hat{\tau}_{wF}= \frac{1}{2} Re \frac{\hat{q}_F |\hat{q}_F|}{\hat{h}^2} C_f,} \end{equation}
having used the reference quantities in table 1. Comparing (4.1) and (4.2), the skin friction coefficient in laminar conditions is 
$C_f=6/Re_F$.
To extend the model to turbulent films, while allowing for recovering (4.1) in laminar films, two constraints should be considered. The first concerns the sign of 
${\tau }_{wF}$, which is given by 
$\hat {q}_F$ in laminar conditions. It is worth noting that in the run-back flow region, as later discussed in the example in figure 16, this quantity is usually negative (i.e. directed upward) since the dominant role of the (positive) shear stress term yields
\begin{equation} {\hat{q}_F=\hat{q}-\tfrac{1}{2} \hat{\tau}_g \hat{h}^2 +\hat{h}<0 \quad \mbox{since} \ \tfrac{1}{2} \hat{\tau}_g \hat{h}^2>\hat{q}+\hat{h}}. \end{equation}
In order to avoid discontinuities in the shear stress in the transition regions, we postulate that the sign of 
$\hat {\tau }_{wF}$ remains dictated by 
$\hat {q}_F$ also in turbulent conditions. The second constraint is to ensure a smooth transition in 
${\tau }_{wF}$ as the flow passes a certain critical Reynolds number. This can be easily ensured if both the mean velocity in the definition of 
$C_f$ and the correlations for 
$C_f$ are solely functions of 
$\hat {q}_F$.
The TTBL proposed in this work is based on the simplest modelling solution respecting these two constraints. First, we keep (4.2) also in turbulent conditions. Second, we adapt the skin friction in the presence of turbulence to a correlation of the form 
$C_f\approx aRe_F^{b}$, with 
$b\approx -1/4$, as encountered in seminal works on turbulent lubrication (Elrod & Ng Reference Elrod and Ng1967; Hirs Reference Hirs1973) and turbulent boundary layer theory (Schlichting & Gersten Reference Schlichting and Gersten2000). Assuming that the transition occurs at 
$Re_*$, the correlation allowing for a continuous transition is
\begin{equation} C_f= \begin{cases} {6}/{Re_F}, & Re_F < Re_*,\\ \left(6 Re^{-3/4}_{*}\right)Re_F^{-0.25}, & Re_F>Re_*, \end{cases} \end{equation}
where the critical Reynolds number here is taken as 
$Re_*=100$. Equations (4.2) and (4.4) can be used for computing the shear stress term 
${\rm \Delta} \hat {\tau }$ in (2.4b).
An alternative formulation could consist in defining the friction coefficient from a reference velocity 
$\hat {q}_T/h$, with 
$\hat {q}_T=\hat {q}_F+1/2\hat {\tau }_g\hat {h}^2$, i.e. including both the first two flow rate contributions in (2.8a–c). In this case, ensuring the aforementioned constraints become more challenging and require the introduction of appropriate blending functions.
4.2. Closures for the advection term
As for the laminar case, the closure of the advection term in (2.5a,b) requires some assumptions on the velocity profile within the film. Self-similar profiles such as the one-seventh-power law have been borrowed from boundary layer theory in the liquid film literature and have shown reasonably good agreement in the prediction of film thickness (Alekseenko et al. Reference Alekseenko, Nakoryakov and Pokusaev1994). More sophisticated eddy viscosity models have been used to analyse phenomena such as gas absorption or heat transfer (Mudawwar & El-Masri Reference Mudawwar and El-Masri1986; Riazi Reference Riazi1996) or the impact of the interface shear stress (Geshev Reference Geshev2014).
Turbulence increases the momentum diffusion, flattening the velocity profile with respect to the laminar profile. This has an impact on the advection term, and possibly on the response of the film thickness. To analyse this impact, we here consider the velocity profile of the falling film portion 
$\hat {u}_F$ (see 2.6) as composed of two contributions, i.e.
\begin{align} \hat u_F(\hat x,\hat y,\hat t)= \hat{u}_L+\hat{u}_T= a_L(\hat x,\hat t) \left[\left(\frac{\hat y}{\hat h}\right)-\frac{1}{2}\left( \frac{\hat y}{\hat h}\right)^2\right]+ a_T(\hat x,\hat t)\left[\left(\frac{\hat y}{\hat h}-1 \right)^{n_T}+1\right], \end{align}
where 
$n_T$ is an integer and odd number. By definition, the turbulent contribution, weighed by the coefficient 
$a_T$, satisfies the boundary conditions for the 
$\hat {u}_F$ term and reproduces boundary layers of different thicknesses with a flat profile above it. The combination of the two terms in (4.5) allows for representing various kinds of departure from the parabolic assumption. The coefficients 
$a_L$ and 
$a_T$ are constrained by the wall shear stress (obtained from (4.2) and (4.4)) and the mass conservation:
\begin{equation} \left. \begin{array}{c@{}} \hat{\tau}_{wF}=\partial_{\hat y} \hat u_F|_{\hat{y}=0}\rightarrow \hat{\tau}_{wF}{\hat h}= a_L + n_T a_T,\\ \displaystyle \hat{q}_F=\int_0^{\hat{h}} \hat u_F \, \textrm{d} \hat y \rightarrow {\hat{q}_F}=\dfrac{1}{3} a_L {\hat h} + \left(\dfrac{n_T}{n_T+1}\right)a_T{\hat h}. \end{array}\right\} \end{equation}The solution of the resulting linear system of equations gives
\begin{equation} \left. \begin{array}{c@{}} \displaystyle a_L=\dfrac{\hat h \hat{\tau}_{wF}n_T}{(n_T+1) c_T}- \dfrac{n_T\hat{q}_F}{\hat h c_T},\\ \displaystyle a_T=-\dfrac{\hat h \hat{\tau}_{wF}}{3 c_T}+\dfrac{ \hat{q}_F}{\hat h c_T}, \end{array}\right\} \end{equation}
where 
$c_T=(2n_T-n^2_T)/(3n_T+3)$. The model is closed for a given 
$n_T$. Regardless of the choice of 
$n_T$, this model recovers the IBL formulation in laminar conditions, for which 
$\hat {\tau }_{wF}=3\hat q_F/\hat h^2$, as this leads to 
$a_T=0$.
While the link between film thickness and flow rate requires the complete numerical analysis of the TTBL model, a first estimation of the range of validity of the model can be inferred from a simple physical constraint: the maximum velocity in (4.5) should be reached at the interface and no other extremes should occur within the liquid film:
\begin{equation} \hat{h} \partial_{\hat{y}}\hat{u}_F =a_L\left(1-\frac{\hat y}{\hat h}\right)+a_T n_T \left(\frac{\hat y}{\hat h}-1 \right)^{n_T-1}\geq 0 \quad \forall \, \hat{y}\leq \hat{h}. \end{equation}
Solving the inequality for any given set of coefficients 
$a_L, a_T$ allows for computing the range of validity of the model for a given 
$n_T$. This range is shown in figure 3(a) for 
$n_T=21$. For each combination of coefficients, it is possible to compute the shape factor of the corresponding profile, defined as
\begin{equation} \varUpsilon=\frac{\hat{h}}{\hat{q}}\int_0^{\hat{h}}\hat{u}^2\, \textrm{d} \hat{y}. \end{equation}
This parameter ranges from 
$1.2$ in laminar conditions to 
$\approx 1$ in case of a flat velocity profile. Using (4.4) and (4.7), it is then possible to analyse how the shape factor changes as a function of 
$Re_F$ for a given 
$n_T$, as well as the maximum 
$Re_F$ tolerated by the model within its range of validity. This is shown in figure 3(b) for 
$n_T=\{7,15,21\}$.
Figure 3. (a) Range of admissible pairs (
$a_L$, 
$a_T$) for a velocity profile in (4.5). (b) Shape factor (4.9) of the velocity profile in (4.5) as a function of the local Reynolds number 
$Re_F$ for 
$n_T=7,15,21$. For each coefficient, the maximum admissible 
$Re_F$ is reported. Plots (c–e) show examples of normalized velocity profiles for the three models at the 
$Re_F$ indicated in each plot, together with the parabolic profile assumed in the IBL model.
Figure 3(c) compares the velocity profiles (normalized to have unitary mean) for 
$n_T=\{7,15,21\}$ with the parabolic assumption at 
$Re_F=354$, i.e. the maximum admissible value for 
$n_T=7$. The same comparison is shown in figure 3(d) at 
$Re_F=930$, the maximum admissible value for 
$n_T=15$, and in figure 3(e) at 
$Re_F=1414$, the maximum admissible value for 
$n_T=21$.
For the purposes of this work, we limit our analysis of the impact of 
$n_T$ on the liquid film modelling to the definition of the upper limit within which such a model is valid. In what follows, a value of 
$n_T=21$ is considered. At the maximum Reynolds number, the resulting velocity profile features a boundary layer of approximately 
$\hat {h}/5$. Compared with the measurements in turbulent falling films (e.g. Mudawar & Houpt Reference Mudawar and Houpt1993), this estimation appears extreme, but well in line with the primary purpose of testing the impact of high turbulence in some portions of the falling film. Nevertheless, it is worth noting that such an extreme is not reached in any of the investigated test cases, among which the highest falling film Reynolds number, in the run-back flow region, is 
$Re_F\approx 600$.
The advection term for the TTBL model, considering 
$n_T=21$, is computed by introducing (4.5) in (2.6) and in the definitions (2.5a,b):
\begin{align} \hat{\mathcal{F}}&=\int^{\hat{h}}_{0} \hat{u}^{2}\, \textrm{d} \hat y =\frac{1}{3} \hat h^3 \hat{\tau}^2_g +\frac{252}{253} a_T \hat{h}^2 \hat{\tau}_g+\frac{5}{12} a_L\hat{h}^2\hat{\tau}_g -\hat{h}^2\hat{\tau}_g+\frac{441}{473} a^2_T \hat{h}\nonumber\\ &\quad + \frac{175}{264} a_L a_T \hat{h}-\frac{21}{11} a_T \hat{h}+\frac{2}{15} a^2_L\hat{h}-\frac{2}{3}a_L \hat{h}+\hat{h}. \end{align}5. The wiping actuators
Following classical modelling strategies of the jet wiping process, the action of the gas jet is modelled via the pressure gradient and the shear stress produced on the liquid film. These two quantities, referred to as wiping actuators, are modelled via experimental correlations for gas jet impinging on a flat (dry) plate (Beltaos Reference Beltaos1976; Tu & Wood Reference Tu and Wood1996; Elsaadawy et al. Reference Elsaadawy, Hanumanth, Balthazaar, McDermid, Hrymak and Forbes2007a; Gosset Reference Gosset2007), with minor adaptations to account for their time dependency, and under the assumption that the dynamics of the liquid film has no influence on their evolution.
Using the reference scales in table 1, the wiping actuators are of the form
\begin{gather} \partial_{\hat x} \hat p_{g}=\partial_{\hat x} [{P}_g(t) f_p(\tilde{x}(t))/(\rho_l \,g)], \end{gather}
\begin{gather}\hat \tau_{g}={T}_g(t) f_\tau(\tilde{x}(t))/(\sqrt{\mu_l\,\rho_l\,g\,U_p}), \end{gather}
where 
$\tilde {x}$ denotes a time-dependent axis accounting for the possible oscillation of the jet, as described at the end of this section. The functions 
$f_p$ and 
$f_\tau$ have range 
$\in [-1,1]$ so that the maximum values from these quantities are defined by the scalars 
${P}_g(t)$ and 
$T_g(t)$. Following the empirical correlation by Tu & Wood (Reference Tu and Wood1996), the pressure distribution for a gas jet impinging on a flat wall is
\begin{equation} f_p(\xi)= \exp(-0.693 \xi^2)+\frac{0.01895 |\xi|}{1+(\xi-1.67489)^2}, \end{equation}
where 
$\xi =x/b$ is a dimensionless coordinate, with the parameter 
$b$ controlling the spreading of the distribution. A qualitative plot of the pressure distribution is shown in figure 1 on the right, with a red dashed line.
As proposed by Beltaos (Reference Beltaos1976), for a stand-off distance 
$Z/d>5$, this parameter can be computed as 
$b=0.125\,Z$. The maximum pressure 
${P}_g$, for a statistically stationary impinging jet, can be computed as 
$P_g=6.5 {P_d d}/{Z}$, where 
$P_d=C_d\,{\rm \Delta} P_N$ is the dynamic pressure at the nozzle outlet and 
$C_d$ is the discharge coefficient taking into account the losses due to friction and separation phenomena in the nozzle chamber. This parameter depends on the nozzle design and is taken as 
$C_d=0.8$ in this work. Considering the reference quantities in table 1, the role of the pressure gradient in the wiping capabilities of the jet is well described by the dimensionless group
\begin{equation} \varPi_g\equiv\frac{P_d d}{\rho_l g Z^2}=C_d \frac{{\rm \Delta} P_N d}{\rho_l g Z^2}, \end{equation}
hereinafter referred to as the wiping number. This number compares the maximum pressure gradient produced by the gas 
$(\sim C_d {\rm \Delta} P_N d/Z^2)$ with the reference (hydrostatic) pressure gradient in the liquid film (
$\rho _l g$).
The distribution of shear stress at the gas–liquid interface is computed following the numerical correlation proposed by Elsaadawy et al. (Reference Elsaadawy, Hanumanth, Balthazaar, McDermid, Hrymak and Forbes2007b). For 
$\xi \geq 0$, this reads as
\begin{equation} f_{\tau}= \begin{cases} {erf}(0.41\xi)+0.54\xi \exp(-0.22\xi^3), & \xi \leq1.73,\\ 1.115-0.24 \ln(\xi), & \xi >1.73. \end{cases} \end{equation}
For 
$\xi <0$, this distribution is mirrored such that 
$f_{\tau }(\xi )=-f_{\tau }(-\xi )$. A qualitative plot of the shear stress distribution is shown in figure 1, with a dash–dotted blue line.
For a sufficiently high Reynolds number in the jet flow, like those considered in this work, the maximum shear stress is computed as 
${T}_g=C_\tau P_d \,d/Z$, with 
$C_\tau =0.067$ (Tu & Wood Reference Tu and Wood1996). Considering the reference shear stress in table 1, the role of the shear stress in the wiping capabilities of the jet is measured by the dimensionless group
\begin{equation} \mathcal{T}_g \equiv \frac{T_g}{Z \sqrt{\rho_l g \mu_l U_p}}= C_d C_\tau\frac{{\rm \Delta} P_N d}{Z \sqrt{\rho_l g \mu_l U_p}}, \end{equation}
hereinafter referred to as the shear number. This number compares the maximum shear stress produced by the gas flow (
$\sim C_{\tau } P_d d/Z$) with the reference shear stress in the liquid film (
$\mu _l U_p/[h]$).
It is worth noting that in the case of wiping of very viscous liquids (such as, e.g. mineral oils or paint) the importance of this number decreases considerably. More information on the scaling laws of the jet wiping process and the typical operating conditions encountered in various industrial processes is presented in Gosset et al. (Reference Gosset, Mendez and Buchlin2019).
The use of correlations for gas jet impinging on dry surfaces has been validated in various studies (see Lacanette et al. Reference Lacanette, Gosset, Vincent, Buchlin and Arquis2006). While the validity of such simplification in time-dependent conditions is certainly questionable, it is important to recall that the modelling of the shear stress produced by an impinging jet on a dry surface is still subject to extensive investigation, and large discrepancies exist in the correlation proposed by various authors. For more details, the reader is referred to Ritcey, McDermid & Ziada (Reference Ritcey, McDermid and Ziada2017).
Finally, concerning the time dependency of the actuators, two possibilities are considered in this work: pulsations and oscillations. In the case of pulsations, the amplitude of the actuators 
$P_g(t)$ and 
$T_g(t)$ are set as harmonics with a mean value equal to the correlations in steady-state conditions and an amplitude of 
$30\,\%$. In the case of oscillations, the amplitudes are left stationary and equal to the correlations previously proposed, while the streamwise variable is taken as 
$\tilde {x}(t)=x-Z\,\tan (W(\theta (t)))$, where 
$\theta (t)$ is the angle of the oscillation with respect to the horizontal, taken as 
$\theta >0$ for a jet deflected upstream (on the run-back flow side), and 
$W(\theta (t))$ is a possible waveform of the oscillation.
Three waveforms are considered. The first is a harmonic oscillation 
$W(\theta (t))=\theta _A\sin (2\pi f t)$, with 
$f$ the perturbation frequency. The others are non-harmonic oscillations biased upstream or downstream. These are constructed by smoothing a square wave signal, which is symmetric around the mean but has a different duration of the positive/negative cycles. In the investigated test cases, a biased oscillation spends 
$80\,\%$ of its period upward or downward. The spatio-temporal evolution of the pressure gradient for an example of each of these oscillatory modes is shown in figure 4. Figures 4(a) and 4(b) show a harmonic and an upward-biased oscillation, while figure 4(c) shows a pulsating perturbation.
Figure 4. Example of two oscillatory waveforms for the jet perturbation: (a) harmonic oscillation; (b) upward-biased (
$\theta <0$) oscillation. The oscillation biased downward (
$\theta >0$) is simply flipped along the 
$\hat x$ axis. The third case (c) is that of a jet pulsation.
These time-dependent jet perturbations were designed to mimic different oscillatory modes in the impinging jet flow. Jet flow oscillations in the wiping process are experimentally investigated in Mendez et al. (Reference Mendez, Gosset and Buchlin2019), while a similar oscillatory mechanism was analysed in Mendez, Scelzo & Buchlin (Reference Mendez, Scelzo and Buchlin2018) on a deformable interface reproducing both fixed and moving surfaces. While it is now known that jet oscillations are coupled to the interface instability, this work aimed at analysing the response of the liquid film to possible jet disturbances, such as oscillations and pulsations, disregarding coupling effects.
6. Numerical methods
6.1. One-dimensional solver for integral models
We here introduce the numerical methods to solve the set of equations (2.4). This is a system of hyperbolic partial differential equations (PDEs) that can be written in the general form
\begin{equation} \partial_t \boldsymbol{V}(x,t)+\partial_{x} \boldsymbol{F}(x,\boldsymbol{V})=\boldsymbol{S}(x,t,\boldsymbol{V}),\end{equation}
with 
$\boldsymbol {V}$ the state vector of the problem, 
$\boldsymbol {F}$ the conservative flux and 
$\boldsymbol {S}$ the source term. In both the IBL and WIBL models, the state vector is 
$\boldsymbol {V}=[h , q]^T$, the flux term is 
$\boldsymbol {F}=[q , \mathcal {F} \beta ]$ and the source term is
\begin{equation} \boldsymbol{S}=\begin{bmatrix}0 \\ (\hat h+ \hat h \partial_{ \hat x \hat x\hat x} \hat h- \hat h \partial_{\hat x}\hat p_g)+{\rm \Delta} {\hat \tau} \end{bmatrix} \frac{\beta}{\varepsilon Re}, \end{equation}
where the coefficient 
$\beta =6/5$ is introduced only for the WIBL and is 
$\beta =1$ otherwise. This coefficient is introduced to account for the contribution 
$\partial _t \hat q/5$ which appears in the shear stress term 
${\rm \Delta} \tau$ (see (B13) in appendix B).
The system of PDEs in (6.1) has been extensively treated in the literature of non-homogeneous shallow water (SW) equations (in which the source term typically accounts for bed topography) and a wide range of suitable finite volume (FV) schemes for their numerical analysis, as described in various textbooks (Toro Reference Toro2001; LeVeque Reference LeVeque2002). Among these, two major classes can be distinguished in the literature: methods based on the (approximated) solutions of the Riemann problem (arising from Godunov's scheme) and methods based on centred fluxes (arising from the Lax–Friedrich scheme). The first class of methods is better suited to handle strong gradients, such as hydraulic jumps, while the second has the advantage of a much lower computational cost (Kurganov & Liu Reference Kurganov and Liu2012; Hernandez-Duenas & Beljadid Reference Hernandez-Duenas and Beljadid2016). Because the investigated simulations do not produce shocks within the space and time domain of interest, this work focused on the second class of methods.
Centred schemes are usually used with a certain amount of artificial viscosity (see Mattsson & Rider Reference Mattsson and Rider2014; Ginting & Mundani Reference Ginting and Mundani2018 and references therein), which can be introduced by a suitable combination of low-order schemes and high-order schemes. This combination is achieved using flux limiters (LeVeque Reference LeVeque2002) to blend a high-order scheme (e.g. Lax Wendroff) in regions where the solution is sufficiently smooth with a low-order scheme (e.g. Upwind or Lax Friedrich) in regions of strong gradients. This approach combines the advantages of the two options: first-order schemes prevent numerical oscillations (dispersion) at the cost of excessively smoothing the solution, while the reverse is true for high-order methods.
A standard FV formulation using explicit methods with three-point stencil in conservative form discretizes (6.1) as
\begin{equation} {\boldsymbol V}_{i}^{k+1} = {\boldsymbol V}_{i}^{k}-\frac{{\rm \Delta} t}{{\rm \Delta} x}[ \boldsymbol{F}^{+}-\boldsymbol{F}^{-}]+ {\rm \Delta} t \boldsymbol{S}^{k}_{i}, \end{equation}
where 
$\boldsymbol {F}^{+}=\boldsymbol {F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i+1}^{k})$ and 
$\boldsymbol {F}^{-}=\boldsymbol {F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i-1}^{k})$ are the fluxes on the right and the left boundaries of each cell. In a flux limiting scheme, these are
\begin{equation} \boldsymbol{F}_i=\boldsymbol{F}^H_i+(\boldsymbol{F}^L_i-\boldsymbol{F}^H_i)\phi_i, \end{equation}
where 
$\phi _i$ is the flux limiting function, 
$\boldsymbol {F}^H$ is the flux calculated from a high-order scheme and 
$\boldsymbol {F}^L$ is the flux calculated from a low-order scheme.
In this work we select the two-step Lax–Friedrich scheme (LxF in Shampine Reference Shampine2005b) as a low-order flux 
$\boldsymbol {F}^L$ and Richtmyer's two-step variant of the Lax–Wendroff method as high-order flux 
$\boldsymbol {F}^H$. An efficient implementation of both schemes in Matlab is provided by Shampine (Reference Shampine2005a), and this work proposed a minor modification to combine the two. These schemes are described in appendix C.
The chosen limiter function is the classical min-mod:
\begin{equation} \phi_i=\max(0,\min(1,\theta_i)). \end{equation}
Here 
$\theta _i=(h_{i}-h_{i-1})/(h_{i+1}-h_{i})$ is the smoothness parameter based on the liquid film thickness as it is customary also in SW problems (e.g. Zhou et al. Reference Zhou, Causon, Mingham and Ingram2001).
The proposed strategy allows for avoiding the calculation of the Jacobian and its eigenmode decomposition, and is therefore computationally cost effective. On the other hand, the numerical diffusion added by the low-order scheme results in the smoothing of the waves in the liquid film. This smoothing becomes more evident as the waves move away from the wiping point from which they originate.
An example of a mesh independency study is shown in figure 5 for an instantaneous thickness and flow rate profile and four different meshes with 
$n_x=\{1940, 2909, 3879, 4848\}$ mesh points. These are computed by setting the number of mesh points 
$n_P$ within the half-width 
$b$ of the Gaussian pressure distribution at the wall from (5.2), so that 
${\rm \Delta} \hat x=b/([x] n_P)$ in dimensionless form. The simulations shown are computed with 
$n_P=\{20,30,40,50\}$, hence ensuring a reasonable accuracy in the calculation of the pressure gradient.
Figure 5. Grid dependency analysis on the film thickness (left) and the flow rate (right). In the zoomed region far downstream of the wiping point (Zoom 2), the effect of numerical diffusion makes the convergence harder than in the wiping region (see Zoom 1). The simulations are carried out with 
${\rm \Delta} P_N=30$ kPa, 
$\theta _A=30^\circ$, 
$Z=15$ mm, 
$d=1.5$ mm, 
$U_p=3$ m s
$^{-1}$, 
$\nu _g=1.5\times 10^{-5}\ \textrm{m}^2\ \textrm {s}^{-1}$, using zinc as a working fluid.
The effect of numerical diffusion, as the number of mesh points is reduced, is evident in figure 5. However, as the focus of this work is placed on the response of the liquid film within a relatively short distance from the introduced perturbation, this is not considered as a limitation.
In all the simulations of this work, the boundary conditions are set as non-reflecting open boundaries while the initial solution is taken from the simplified one-dimensional formulation. The simulation is run until a fully periodic response is produced in the film before the data for post processing is acquired. The time step is taken by setting 
${\rm \Delta} t=0.4 {\rm \Delta} x$, i.e. assuming a Courant–Friedrichs–Lewy (CFL) number equal to 
$u_w\,{\rm \Delta} t/{\rm \Delta} x=0.8$ for waves travelling at 
$u_w\approx 2$, that is, twice the substrate speed. Such an estimation revealed to be rather conservative.
6.2. Direct numerical simulations and validation of the long-wave formulation
Before considering the jet wiping problem, we analyse the validity of long-wave formulations on a much simpler test case, namely the flow of a wavy liquid film over an upward-moving substrate. This configuration is relevant to describe the dynamics of the final coating much more downstream the wiping region, where the pressure gradient and the shear stress produced by the impinging gas jet vanish.
Although this test case is too simple for complete validation of the models, which is out of the scope of this work, we here focus on the validity of the long-wave formulation of the jet wiping problem at large Reynolds numbers. Moreover, we analyse the validity of the proposed Shkadov-like scaling by considering two liquids with largely different properties operating at the same dimensionless thickness 
$\hat {h}$ and rescaled Reynolds number 
$\delta =\varepsilon Re$. The liquids considered are water and zinc.
The validation has been carried out using direct numerical simulations of the two-phase flow using the VOF method, in which the gas–liquid interface is tracked on a fixed grid. Surface tension is accounted for through the continuum surface force method (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992), in which the surface force due to capillarity is converted into a volume force that acts across the interface thickness. The computational domain is shown in figure 6(a) along with the required boundary conditions. The liquid film has a mean (and initial) thickness 
$h_o$.
Figure 6. (a) Schematic of the flow configuration used for validation purposes using the OpenFoam solver InterFoam: flow domain and zoom on the near-wall mesh. (b) Snapshot of the thickness evolution for the four meshes in table 2.
The computational domain is rectangular, with a dimensionless length 
$L_x=8400 h_o$ for water and 
$L_x=12300 h_o$ for zinc in the streamwise direction, whereas 
$L_y=7.8 h_o$ for water and 
$L_y=7.5 h_o$ for zinc in the cross-stream direction. A perturbation of the film flow rate is introduced at the inlet of the domain via a pulsation of the film velocity (as in Doro & Aidun Reference Doro and Aidun2013). At the liquid inlet (with thickness held constant and equal to 
$h_o$), the film velocity profile is prescribed as
\begin{equation} \hat{u} (\hat{x}=0,\hat y)=[\hat y (\tfrac12 \hat y - \hat h_o )+1][1+q_A\sin(2\pi\hat{f}\hat t)], \end{equation}
where 
$\hat f=f [t]$ is the dimensionless oscillation frequency, 
$\hat t$ is the dimensionless time, and assuming that the 
$x$ axis points out in the same direction as the substrate velocity 
$U_P$. The corresponding flow rate per unit width is therefore
\begin{equation} \hat{q}=[\tfrac13 \hat{h}^3-\hat h][ 1+q_A\sin(2\pi\hat{f}\hat t) ]. \end{equation}For these computations, the interFoam solver of the finite volume code OpenFoam is used. This solver has been extensively validated in the computational fluid dynamics literature (Deshpande, Anumolu & Trujillo Reference Deshpande, Anumolu and Trujillo2012) and in various studies on falling liquid films (e.g. Gao et al. Reference Gao, Morley and Dhir2003; Doro & Aidun Reference Doro and Aidun2013; Dietze et al. Reference Dietze, Rohlfs, Nährich, Kneer and Scheid2014) as well as co-current and counter-current gas–liquid flows (Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2013). The solver assumes incompressible and isothermal flow and features an interfacial compression flux term that activates at the interface to mitigate the effects of numerical smearing of the gas–liquid boundary.
The liquid volume fraction 
$\alpha$ is fixed at the inlet, together with a Neumann condition for pressure. On the inlet boundary located in the gas phase (
$x _0=0$, 
$h_o \leq y \leq 7.8h_o$) and on the right-hand side of the gas boundary (
$y=7.8 h_o$, 
$0 \leq x \leq L_x$), 
$\alpha$ is fixed to 0, with a zero derivative for the velocity and a fixed total pressure. At the wall (
$y=0$, 
$0 \leq x \leq L_x$), a no-slip condition is prescribed (
$\hat u(\hat y=0)=1$), with a zero flux for 
$\alpha$ and a fixed pressure condition. At the outlet, a zero gradient is set for both the velocity and the liquid volume fraction. The flow field is initialized with the nominal thickness 
$h_o$. The velocity profile within the liquid is initially parabolic (
$\hat {t}=0$ in (6.6)) while the velocity is set to 0 in the gas phase. In agreement with (3.14), we consider cases with 
$\hat {h}\ll 1$ such that the waves propagate upstream, which is in the direction of the substrate motion.
The simulations are carried out using a second-order backward Euler scheme in time for the transient term, and second-order discretization schemes for the convective (van Leer scheme for the 
$\alpha$ transport equation), diffusive and pressure terms. The coupling between pressure and velocity is solved using a standard PISO algorithm. The time step was set adaptatively, based on a maximum value of 0.3 for the global CFL number in the 
$\alpha$ equation. This leads to time steps of the order of 
$4.5 \times 10^{-6}$ s with water and 
$2.7 \times 10^{-6}\,{\rm s}$ with zinc.
In order to evaluate the influence of mesh density on the results, simulations are performed on four different grids with an increasing mesh density (see table 2): the streamwise cell size 
${\rm \Delta} x$ is varied between 
$0.11 h_o$ and 
$0.39 h_o$, and the cross-stream cell size 
${\rm \Delta} y$ between 
$0.022$ and 
$0.078 h_o$. For these tests, water was used as the working fluid, and the length of the domain was reduced to save computational time (
$L_x=1500 h_o$). The substrate speed is fixed to 1 m s
$^{-1}$. The results in figure 6(b) show that the selected mesh densities are sufficient to capture the sinusoidal waves that form shortly after the inlet, and that beyond the density of mesh M2, the thickness profiles are almost insensitive to the size of the cells. The meshes used in the validation process have cell sizes of the order of 
${\rm \Delta} x=0.23 h_o$ and 
${\rm \Delta} y=0.046 h_o$, resulting in grids of 3.3 million for the test case with water and 4.6 million for the one with zinc.
Table 2. Mesh densities for the sensitivity study.
7. Results
7.1. Validation test cases
The physical parameters and the operating conditions for the two simplified cases with water and zinc are recalled in table 3. Both test consider a dimensionless liquid thickness of 
$h_o/[h]=0.2$, perturbed at the inlet with a flow rate pulsation with amplitude 
$q_A=0.2$ and dimensionless frequency 
$\hat {f}=0.05$. The substrate velocity is taken as 
$U_p=1$ m s
$^{-1}$ for both cases. These two liquids differ by one order of magnitude in surface tension and almost an order of magnitude in density and dynamic viscosity. However, their kinematic viscosity is comparable, resulting in similar time scales (cf. table 1), and similar rescaled Reynolds number 
$\delta =\varepsilon Re$.
Table 3. Physical parameters and operating conditions of the two test cases used for validation purposes.
An instantaneous thickness profile is shown in figure 7(a) for both cases, comparing the IBL, WIBL and OpenFoam simulations. Since 
$Re_F\approx 95$ in zinc and 
$Re_F\approx 65$ in water, the TTBL recovers the IBL model and its results are not shown. The thickness profiles are perfectly overlapping, demonstrating the validity of the integral formulation and the numerical methods, as well as the capability of the long-wave formulation to model the flow.
Figure 7. (a) Instantaneous thickness profile for the liquid film thickness of the liquid zinc (left) and the water (right) test cases. (b) Comparison of the velocity profiles extracted from the zinc simulations in the OpenFoam (OF) and IBL/WIBL simulations. The profiles are extracted from the wave maxima (location 1) and wave minima (location 2) in (a).
Table 4 collects all results in terms of flow rate and thickness maxima and minima, while figure 7(b) shows the velocity profile underneath a maximum (1) and a minimum (2) of the film waves in the simulations with zinc. The location in which the profiles are recomputed is indicated in figure 7(a). For both IBL and WIBL, these profiles are reconstructed from the results of the simulation (
$\hat {h},\hat {q}$). In the IBL model this is a straightforward implementation of (3.6); in the case of the WIBL model this involves all the relations in table 6 with (3.1) and (2.6). For the WIBL model, reconstructing the velocity profile is an ill-posed problem, since the closure of the model (i.e. the expression for 
${\rm \Delta} \hat {\tau }$ in (B13)) assumes that 
$\delta \ll 1$ while in this case 
$\delta \sim 75$. Nevertheless, in these simple examples, the 
${\rm \Delta} \tau ^{(1)}$ term is small enough to let the WIBL converge on the IBL despite the large 
$\delta$, and the reconstructed profile reflects this perfect convergence. As later discussed in § 7.4, this does not happen for the jet wiping configurations at higher Reynolds numbers, as the higher-order corrections of the velocity profile become more important than the zeroth order.
Table 4. Summary of the results in terms of flow rate and thickness maxima/minima (subscript min/max) and wavelength (
$\hat {\lambda }$) for the IBL/WIBL models and the OpenFoam (OF) simulation.
Finally, it is worth highlighting excellent agreement between the integral models and the direct numerical simulation calculations, which reveals parabolic velocity profiles. This result shows the important effect of the substrate motion as compared with the classic falling film problem, where the departure from the parabolic profile occurs at a much lower Reynolds number than the values considered here (see Denner et al. Reference Denner, Charogiannis, Pradas, Markides, van Wachem and Kalliadasis2018).
7.2. The relative contribution of forces
This section focuses on the relative importance of all the terms in the integral momentum formulation (2.4b) as the wiping strength (as measured by the wiping number 
$\varPi _g$) is increased. The WIBL and TTBL models are considered. As in the previous section, both liquid zinc and water are analysed as working fluids.
The liquid properties are the same as the previous section, as listed in table 3. In the case of zinc, the wiping conditions are taken to be representative of an industrial galvanizing line, with a nozzle having an opening of 
$d=1.5$ mm and stand-off distance 
$Z=15\,\rm {mm}$. The substrate velocities 
$U_p=\{1,2,3\}$ m s
$^{-1}$, corresponding to 
$Re=\{478,1352,2483\}$ and 
$\delta =\{74,263,554\}$, are considered and, for each of these, the pressure in the nozzle is varied in the range 
${\rm \Delta} P_N=[3,40]$ kPa. This leads to wiping numbers in the range 
$\varPi _g=[0.7,2.9]$ and a shear stress number in the range 
$\mathcal {T}_g=[5,70]$.
In the case of water, the nozzle parameters 
$Z$ and 
$d$ are the same as for the zinc cases. The substrate velocity is reduced to 
$U_p=\{0.2,0.3,0.4\}$ m s
$^{-1}$, corresponding to 
$Re=\{28,52,80\}$ 
$\delta =\{4,8.4,14.3\}$, and the nozzle pressure decreased down to the range 
${\rm \Delta} P_N=[0.7,1.5]$ kPa. This leads to much lower wiping numbers, in the range 
$\varPi _g=[0.3,0.8]$, but a comparable shear stress number, in the range 
$\mathcal {T}_g=[30,80]$. The simulated wiping conditions fall within the operational range of the Essor VKI facility (Buchlin Reference Buchlin1997).
In all the investigated configurations, the harmonic oscillations of the jet have an amplitude of 
$\theta _A=10^{\circ}$. It is worth noting that in the Shkadov scaling, this leads to different amplitudes of the oscillations along the streamwise direction (
$\hat {x}$) as the jet stand-off distance is the same while the streamwise length scale is not.
Three dimensionless frequencies are considered, i.e. 
$\hat {f}=\{0,0.055,0.16\}$. The zero frequency simulates the process in steady-state conditions, yet accounting for the role of inertia and surface tension. Several preliminary tests in this configuration confirm that the flow is absolutely stable, meaning that initial disturbances in the film move away from the domain (towards 
$\hat {x}\rightarrow \infty$ if the disturbance is located in the run-back flow; towards 
$\hat {x}\rightarrow -\infty$ otherwise) leaving the steady-state solution unvaried. The highest frequency 
$\hat {f}=0.16$, as further discussed in § 7.3, is damped by the liquid film and results in no appreciable undulation. On the contrary, the frequency 
$\hat {f}=0.055$ is in the range of maximum receptivity of the liquid film and produces the largest waves.
The time-averaged coating thickness downstream of the wiping is shown in figure 8 for both liquids. The left column refers to the wiping conditions in zinc, and the right column to the wiping conditions in water. From top to bottom, the velocity of the substrate is increased, and the corresponding rescaled Reynolds number 
$\delta =\varepsilon Re$ is indicated. Each figure compares the prediction of the ZO model presented in § 3.1. This comparison is made with (black dashed lines) and without (continuous red line) the shear stress in the model. The results of the ZO models are well described by power laws of the form 
$a \varPi ^{b}_g$, reported in the figures.
Figure 8. Time-averaged thickness in the final coating film for three rescaled Reynolds numbers (indicated in each plot) and various oscillation frequencies (
$\hat {f}=\{0,0.055,0.16\}$). The black dashed line and the red continuous line correspond to the prediction from the zero order, the first taking into account the contribution of the shear stress, the second setting the shear stress to zero. Black markers refer to simulations using the WIBL model; white markers refer to simulations using the TTBL model.
In the presence of shear stress, the power correlation changes slightly with the substrate speed, while this remains unaltered if the shear stress is removed. This result is in remarkable agreement with the experimental correlations presented in previous experimental works (Gosset et al. Reference Gosset, Mendez and Buchlin2019; Mendez et al. Reference Mendez, Gosset and Buchlin2019). It is interesting to observe that these experimental works were carried out on a much more viscous mineral oil, producing similar wiping numbers (in the range 
$\varPi _g=[0.1,0.8]$) but much lower shear stress numbers (in the range 
$\mathcal {T}_g=[1,8]$). This highlights the role of the shear stress in the wiping process for liquids with low kinematic viscosity such as zinc or water and confirms the experimental observation that the wiping of highly viscous liquids is mostly governed by the dimensionless group 
$\varPi _g$.
Concerning the role of the substrate speed on the mean film, the ZO model predicts a negligible impact in all the test cases, while the WIBL and TTBL models reveal a discrepancy that becomes more important at lower wiping numbers and large Reynolds numbers. Overall, the WIBL and TTBL models are in good agreement in all the simulations analysed (at small Reynolds numbers, these are indistinguishable). The results from these models at 
$\hat {f}=0$ (trends with square markers) show that surface tension and advection produce a significant departure from the wiping curve obtained by the simplified model at 
$\varPi _g\rightarrow 0$ while the agreement is asymptotically reached at 
$\varPi _g\rightarrow \infty$.
It is worth noting that the cases with jet oscillation yield a higher mean coating thickness regardless of the oscillation frequencies. This result is particularly interesting if one considers that the case at 
$\hat {f}=0.16$ yields no undulation in the final coat as later discussed in § 7.3. Yet, the mean thickness increases, as if the oscillation spreads momentum and results in an effective distribution that is closer to the time-averaged profiles. A similar phenomenon is also observed by Lunz & Howell (Reference Lunz and Howell2018), who studied the response of a liquid film to an oscillatory pressure source and discussed the calculation of the effective pressure.
Focusing on the wiping cases with zinc, figure 9 presents several instantaneous profiles for six representative test cases taken from the results in figure 8. These include the lowest (
$Re=478$) and the largest (
$Re=2483$) Reynolds numbers and the three perturbation frequencies, keeping the wiping number at 
$\varPi _g=1.856$. In each of the six panels, the first plot compares the instantaneous dimensionless film thickness profiles for WIBL and TTBL models. The second plot compares the advection terms 
$\partial _x \mathcal {F}$ profiles for both models while the third collects the contributions to the wall shear stress. The term 
${\rm \Delta} \hat {\tau }^{(0)}$ is the ZO term from the WIBL, which is the one in the IBL model. The term 
${\rm \Delta} \hat {\tau }^{(1)}$ is the first-order contribution that distinguishes the WIBL from the IBL. Finally, the term 
${\rm \Delta} \hat {\tau }_T$ is the wall shear stress term from the TTBL model.
Figure 9. Results from the simulations of the jet wiping in galvanizing conditions using the WIBL and TTBL models for 
$\varPi _g=1.856$, two Reynolds numbers (
$Re=478$ on the left column and 
$Re=2483$ on the right) and three oscillation frequencies (
$\hat {f}=\{0,0.055,0.16\})$ along the rows. In each panel the first plot from the left shows an instantaneous film thickness, the second and the third the corresponding advection and shear stress terms, respectively.
As expected, no differences are observed between the three models at the lowest 
$Re$, except for the cases at the highest perturbation frequency 
$\hat {f}=0.16$. At 
$Re=478$, the wall shear stress terms have a larger contribution to the liquid film dynamics, especially in the proximity of the wiping region. In this condition, the first-order term 
${\rm \Delta} \tau ^{(1)}$ appears to have a negligible contribution; hence, the WIBL model corresponds almost everywhere to the IBL model (not shown). In the case at 
$Re=2483$, the contribution of this term increases but remains less important than the advection term, which mostly dominates the dynamics of the liquid waves. The TTBL model departs from the laminar models in the run-back flow region, while no appreciable difference is observed in the final coating region as this is characterized by 
$Re_F<100$. This further highlights the convective nature of the problem with two opposite characteristic lines: the dynamics in the final coating film appears to be insensitive to the dynamics of the run-back flow region. Finally, in terms of the frequency response of the liquid coat, the three models reveal that the perturbation frequency of 
$\hat {f}=0.16$ is too high to generate any appreciable wave in the final coat. The influence of the modelling strategy on the harmonic response of the flow is discussed in § 7.4; the next section focuses on the harmonic response of the film considering only the WIBL model.
7.3. The frequency response of the liquid film
This section focuses on the frequency response of the liquid film subject to different kinds of jet perturbation, producing pressure gradient evolutions of the form described in figure 4. For a wiping number 
$\varPi _g=1.2$ and rescaled Reynolds number 
$\delta =554$, figure 10 shows several contour maps of the mean-centred thickness 
$\tilde {h}=\hat {h}(\hat {x},\hat {t})-\bar {h}(\hat {x})$, obtained by subtracting the temporal average 
$\bar {h}(\hat {x})=({1}/{T})\int ^T_0 \hat {h}(\hat {x},\hat {t})\textrm {d}\hat t$, with 
$T$ the dominant wave period. Four jet perturbations are considered, namely three oscillations (harmonic, upward biased and downward biased) and a harmonic jet pulsation (cf. figure 4). For each of these, the frequencies considered are 
$\hat {f}=\{0.02,0.05,0.08\}$. For the jet oscillation test cases, a white dashed line indicates the evolution of the impingement point, i.e. the region of maximum pressure and zero gas shear stress at the interface.
Figure 10. Contour maps of the coating thickness 
$\tilde {h}(\hat {x},\hat {t})$, zero-mean shifted in time, as a function of dimensionless space and time. Four kinds of jet perturbation are considered: harmonic oscillations (first row), upward-biased (second row) and downward-biased (third row) oscillations and jet pulsations (fourth row). The perturbation frequencies are taken as 
$\hat {f}=0.02$ (first column), 
$\hat {f}=0.05$ (second column) and 
$\hat {f}=0.08$ (third column). All cases refer to galvanizing conditions with 
$\varPi _g=1.2$ and 
$\delta =554$.
The liquid film response to harmonic oscillation (first row of figure 10) is discussed first. At 
$\hat {f}=0.02$, the coating thickness is characterized by wave peaks of 
$\tilde {h}\approx 0.1$. The characteristic lines tracing the propagation of the waves clearly show that these originate below the average impact point, at 
$\hat {x}\approx -1$, that is, in the region normally belonging to the run-back flow. Within the range 
$\hat {x}\in [-2,2]$, spanned by the jet during the oscillation, regions of liquid film accumulation (
$\tilde {h}>0$) and depletion (
$\tilde {h}<0$) alternate harmonically as the coating film follows the jet oscillation. Within the region intersected by the jet oscillation, the propagation speed of the coating waves is not constant and strongly influenced by the displacement of pressure gradient and shear stress. Outside this region, the wave propagation speed remains constant and approximately equal to the substrate speed for 
$\hat {x}>2$. The liquid meniscus follows the displacement of the wiping region and the contour map of the mean-shifted thickness 
$\tilde {h}$ is almost symmetric along 
$\hat {x}=0$.
This is not the case at higher frequencies (
$\hat {f}=0.05$ and 
$\hat {f}=0.08$), in which the disturbances in the run-back flow appear comparatively much lower and the waveform of the liquid changes significantly. To better analyse the difference in the wave formation mechanism, we consider the response of the liquid film thickness 
$\hat {h}^*(\hat {t})$ in the impact point 
$\hat {x}^*(\hat {t})$, defined as the point in which the gas pressure is maximum. Figures 11(a) and 11(b) show, respectively, the evolution of 
$\hat {h}^*(\hat {t})$ as a function of the dimensionless time, scaled by the wave period, and its phase portrait with 
$\hat {x}^*(\hat {t})$, describing the oscillation of the impact point. For plotting purposes, the curves in figure 11(a) are shifted to have matching wave peaks, while the phase portraits are constructed by mean shifting both signals and normalizing with respect to their peak-to-peak amplitude.
Figure 11. Plots (a) and (c) show the film thickness at the impact point 
$\hat {h}^*(\hat {t})$ as a function of the period normalized time for three dimensionless frequencies, considering jet oscillations (a) and jet pulsations (c). Plots (b) and (d) show the phase portrait linking 
$\hat {h}^*(\hat {t})$ to the jet disturbance: this is the time varying location of the impact 
$\hat {x}^*(\hat {t})$ in case of an oscillation (in b) and the time varying maxima 
$\max(\partial _{\hat {x}\hat {p}_g})$ in case of a pulsation (in d).
At the frequency of 
$\hat {f}=0.02$, the film thickness at the impact point remains almost sinusoidal and has a constant phase delay of approximately 
$\pi /4$. This phase delay is due to the response time of the liquid film as the jet moves towards the run-back flow, encounters a region of higher thickness, and imposes a thickness reduction. At 
$\hat {f}=0.05$ and 
$\hat {f}=0.08$, the response 
$h^*(t)$ is no longer sinusoidal and the film does not have enough time to allow for the wiping: during the ascending phase of the jet oscillation (denoted as A in figure 11a), a portion of un-wiped liquid is dragged from the run-back flow region and pushed towards the final coating region, from which it continues its upwards evolution at the speed of the substrate. As a result, the wave peak reached in the case of 
$\hat {f}=0.05$ is about 20 times higher than in the case of 
$\hat {f}=0.02$. These results confirm the existence of the mechanisms for the wave formation originally postulated in the previous experimental investigation by the authors (Mendez et al. Reference Mendez, Gosset and Buchlin2019) and which hereinafter is referred to as mechanism A.
The non-harmonic cases in the second and third rows of figure 10 show that both the shape and the maximum amplitude of the resulting waves are strongly influenced by the waveform of the jet oscillation. Nevertheless, the key observation from the harmonic case applies: the characteristic lines through which the waves evolve (both towards the run-back flow and towards the final coating film) originate below the mean wiping point.
It is now instructive to consider the pulsating test cases in the last row of figure 10. In this case, as the impingement point is fixed in time, the waves originate at the wiping line 
$\hat {x}=0$ and move at a constant speed in both downward and upward directions. Figures 11(c) and 11(d) show, respectively, the time evolution of the thickness at the impact point and its phase portrait with the maximum pressure gradient during the jet pulsation. The same plotting adjustments of figure 11(a,b) in terms of shifting and normalization are adopted. Regardless of the pulsation frequency considered, the thickness at the impinging point remains overall sinusoidal, with the phase delay converging towards 
$\pi /2$, which is a perfect quadrature. This second mechanism of wave formation is herein referred to as mechanism B.
To conclude this section, figure 12 collects the frequency response of the film coating in both the final coating region and the run-back flow region for the four perturbations considered and for four combinations of wiping number 
$\varPi _g$ and rescaled Reynolds number 
$\delta$. All the test cases refer to galvanizing conditions, and the selected pairs 
$(\varPi _g$,
$\delta )$ are indicated in the legend. For each of the wiping conditions, the plot shows the dependency of the wave amplitude, measured in terms of the standard deviation 
$h_{\sigma }$ to average 
$\bar {h}$ ratio, over the dimensionless frequency. For the final coating film, these curves are computed at a location 
$\hat {x}=-15$ while 
$\hat {x}=15$ is considered for the run-back flow.
Figure 12. Amplitude of the coating waves in terms of 
$h_{\sigma }/\bar {h}$ as a function of the dimensionless perturbation frequency for four wiping conditions, indicated in the legend, and for the four jet perturbations previously considered. The curves on the left column are computed in the final coating region (at 
$\hat {x}=-15$) while the curves on the right are computed in the run-back flow region (at 
$\hat {x}=15$). All the tests refer to galvanizing conditions.
As expected from the previous analysis of the contour maps, the amplitude of the coating waves in the case of an oscillating jet (mechanism A) is significantly larger than in the case of a pulsating jet (mechanism B). Moreover, while mechanism A shows a region of strong receptivity in the range of dimensionless frequencies 
$\hat {f}=[0.03{-}0.08]$, the liquid film behaves as a low pass filter with respect to mechanism B. Regardless of the mechanism, no coating waves can be expected for perturbations at 
$\hat {f}>0.2$.
By comparing the modulation curves for the two regions of the coating flow, the lower portion of the receptivity band, (say 
$\hat {f}\approx [0.03{-}0.05]$) appears of significant interest. This range of frequency is also strongly present in the run-back flow, while the second portion (say 
$\hat {f}\approx [0.05\text{--}0.08]$) is more attenuated. In the coupling mechanism described in previous studies (Pfeiler et al. Reference Pfeiler, Eßl, Reiss, Riener, Angeli and Kharicha2017; Mendez et al. Reference Mendez, Scheid and Buchlin2017b, Reference Mendez, Gosset and Buchlin2019), the jet oscillation is sustained by the waves in the run-back coating flow. These waves appear to be the major cause of the unstable interaction between the two flows.
Concerning the range of possible frequencies, it is interesting to report that 
$\hat {f}=0.05$ – at which maximum amplitude of the coating waves can be expected – corresponds to wavelengths of the order of 
$\lambda =25\text {--}35$ mm depending on the wiping conditions. This range is in agreement with undulation defects observed in several galvanizing lines (Pfeiler et al. Reference Pfeiler, Eßl, Reiss, Riener, Angeli and Kharicha2017). Moreover, the results show that the undulation amplitude produced by an oscillating jet (mechanism A) increases at larger substrate speeds – a fact also in line with industrial observations – while it remains rather insensitive for a pulsating jet (mechanisms B).
As to the role of the wiping number, its impact on the coating wave amplitude is a more difficult interpretation, especially if one considers that the wiping of liquids with low kinematic viscosity such as liquid zinc or water is also strongly influenced by the shear stress, as shown in figures 8. While increasing the wiping number results in a thicker run-back flow, the main role of the shear stress is that of smoothing the wiping meniscus and, thus, the transition from the final coating film to the run-back flow. The smoother this meniscus is, the smaller the thickness gradient encountered by the jet during oscillation, and, hence, the lower the impact of mechanism A. These results also show that the interaction between the gas and the liquid film is strongly influenced by the waveform of a possible jet oscillation: while the receptivity range of the film remains unaffected, the amplitude of the undulation is significantly larger in a harmonic oscillation than in non-harmonic oscillations.
Finally, it is interesting to compare the dimensionless modulation function obtained in the galvanizing conditions with those obtained for the wiping of water. Only the response to harmonic oscillation is considered, and shown in figure 13 for four wiping conditions. Despite the largely different dimensionless numbers controlling the process, the maximum undulation amplitude in the final coating film is produced in a similar range of dimensionless frequencies 
$\hat {f}=[0.03{-}0.08]$ – i.e. the band dominated by mechanism A. Given the largely different properties of the two liquids, these results show that the Shkadov-like scaling well describes the physical phenomena governing the maximum receptivity of the liquid film.
Figure 13. Same as in figure 12 but considering the wiping of water instead of liquid zinc, and harmonic oscillations only.
7.4. The influence of the modelling strategy
It is finally of interest to analyse the impact of the modelling strategy on the results previously obtained by comparing the IBL, WIBL and TTBL models. We here consider the response of the liquid film to harmonic jet oscillations in galvanizing conditions. Figure 14 compares the frequency response from the three models on the two extreme cases from the first row of figure 12: these are 
$\varPi _g=1.2$ and 
$\delta =554$ (largest response, filled markers) and 
$\varPi _g=2.3$, 
$\delta =74$ (smallest response, empty markers). The impact of the model becomes more pronounced as larger waves are considered, with the TTBL model predicting significantly larger waves both in the final coating film and in the run-back flow. It is nevertheless interesting to observe that the range of maximum receptivity of the liquid film remains unvaried, with the three curves showing the same qualitative behaviour.
Figure 14. Same plot as the first row of figure 12, considering the two extreme cases 
$\varPi _g=1.2$, 
$\delta =554$ (largest response, filled markers) and 
$\varPi _g=2.3$, 
$\delta =74$ (smallest response, empty markers) comparing the results from the IBL, WIBL and TTBL models.
Figure 15 shows five snapshots of the film thickness computed by the three models, together with the instantaneous pressure gradient profile, for two test cases with oscillation frequency 
$\hat {f}=0.04$ (on the top) and 
$\hat {f}=0.06$ (on the bottom) and wiping at 
$\varPi _g=1.2$ and 
$\delta =553$. These yield the largest waves in the final coating film mostly originated by mechanism A described in the previous section.
Figure 15. Sequence of snapshots for test cases with 
$\varPi _g=1.2$ and 
$\delta =554$ and an harmonic oscillation of the impinging jet at 
$\hat {f}=0.04$ (top) and 
$\hat {f}=0.06$ (bottom). Each snapshot plots the film thickness for IBL (dashed black line), WIBL (continuous blue line) and TTBL (dotted red line), together with the pressure gradient imposed by the impinging jet (light blue dotted line). The thickness axis is on the top of each plot; the pressure gradient axis is on the bottom. The time step is indicated in each plot. Animations of both cases are provided as supplementary movies 1 and 2, available at https://doi.org/10.1017/jfm.2020.1075.
The sequence of five snapshots captures one period of this mechanism, starting from its most downward position (first column). In the second column, the effect of the wiping is visible while the snapshot in the third column captures the formation of the wave in the final coat, as a local minimum of thickness is dragged upward. In the snapshots of the fourth and fifth columns, the jet is in its descending phase, impinging on a much thicker film, and mechanism A begins its next period.
It is interesting to observe that these waves become strongly asymmetric soon after their formation, with a steeper gradient on their tail. This asymmetry changes as the wave evolves downstream under the action of the shear stress, which imparts higher advection velocity to regions of higher thickness. While the differences between the laminar models are minor, the turbulence model leads to largely different shapes of the waves in the final coating film. Since in this region the TTBL model recovers the IBL, this difference is linked to the discrepancy in the prediction of the thickness of the run-back flow, where the TTBL yields a much thicker film.
Focusing on the influence of the turbulent modelling, figure 16 further analyses the test case with 
$\hat {f}=0.04$ in the upper row of figure 15. The first plot on the left shows the maxima, minima and mean distribution of the term 
$\hat {q}_F Re$, the absolute value of which controls the transition to turbulence. This term is negligible from 
$\hat {x}>15$, suggesting that a laminar model is appropriate to simulate the final coating flow. On the other hand, this term is large in the run-back flow and, most importantly, negative: as described by (4.3), this is a consequence of the large and positive contribution of the shear stress (
$1/2 \hat {h}^2\hat {\tau }_g Re$), depicted in the figure on the right. Mechanism A for the wave formation propagates the discrepancy in the run-back flow towards the final coat.
Figure 16. The plots on the left shows the minima, mean and maximum distribution for the flow rate contributions 
$\hat {q}_F$ and 
$\frac {1}{2}\hat {h}^2\hat {\tau }_g$, multiplied by the global Reynolds number 
$Re$: the absolute value of these quantities leads to the Reynolds numbers in (2.8a–c). The plots on the right show the velocity profiles from the IBL and TTBL under wave peaks in the locations 
$\hat {x}=-50,-15,15$.
Finally, the plots on the right of figure 16 show three instantaneous velocity profiles extracted on a wave peak located at 
$\hat {x}=\{-50,-15,15\}$ for both the IBL and TTBL models. Far from the impact point, in the final coating, the velocity profile appears almost linear due to the thin thickness of the layer. Close to the impact point, in the run-back flow region, the velocity profile appears almost linear due to the strong influence of the gas shear stress, although the discrepancy between the two models increases. Finally, it is worth reporting that in this condition, the rescaled Reynolds number is too high to allow for reconstructing the velocity profile in the WIBL model as its first-order corrections, designed for 
$\delta \ll 1$, become non-physical. In particular, at large 
$\delta$, the coefficients 
$a_j$ with 
$j>0$ in table 6 lead to extremely large contributions to the velocity profile. It appears thus surprising that the predicted film response is still closely matching the IBL, as the inconsistency in the large 
$a_j$'s is somewhat mitigated by the shear stress term in (B12). To conclude, the mechanisms for wave formation revealed in this work, and their range of maximum receptivity are qualitatively independent of the models implemented.
8. Conclusions
We have presented an extension of classical low-dimensional models for falling liquid films to the jet wiping problem, tailoring accordingly the Shkadov-like scaling. This process consists of using an impinging gas jet to control the thickness of a liquid film on a moving substrate and is characterized by an unstable interaction between the gas and the liquid film that limits the achievable coating uniformity.
The investigated integral models allow for simulating this complex interaction in industrially relevant conditions with minor computational costs and, thus, enable insights on the process dynamics that would otherwise not be possible using high fidelity simulations. The proposed models extend the modelling strategies commonly used in falling liquid films to a more complex scenario that includes the motion of the substrate and the presence of an imposed pressure gradient and shear stress distribution – in this work simulating an unstable impinging jet. The extended models, the self-similar IBL and WIBL models, were extensively described and framed along with the classic ZO formulation encountered in the literature of the jet wiping process. Moreover, an extension of the IBL model, referred to as the TTBL model, has been proposed to account for the impact of turbulence in the liquid film response.
The numerical implementation of these models has been successfully validated via direct numerical simulations using the VOF method in OpenFoam, considering the simplified test case of a pulsating liquid film evolving along with a moving interface. These models were then used to study the response of the liquid film to various kinds of perturbation in the jet flow, including harmonic and non-harmonic oscillations and pulsations. The analysis of the relative influence of all the terms in the equations reveals that the nonlinear advection term dominates over a wide range of wiping conditions and frequencies of the perturbation. Moreover, by analysing the wiping process in galvanizing conditions (using zinc) and in laboratory conditions (using water), it is shown that the Shkadov-like scaling reveals an interesting similarity of the frequency response of the liquid coat. In the simplest case of a film flowing over an upward-moving surface, the similarity between the two configurations applies to the entire film evolution.
Two main mechanisms for the formation of waves in the coating film downstream of the wiping region were identified. The first mechanism, referred to as mechanism A, is inherently linked to the presence of a wiping meniscus and to the fact that, during its oscillation, the jet drives liquid upwards from the thicker region below the wiping point. The amplitude and shape of the coating waves produced by this mechanism were shown to be linked to the waveform of the jet oscillation. The second mechanism, referred to as mechanism B, is related to the local variation of the pressure gradient and shear stress, as a result of unsteadiness in the jet flow. The dimensionless transfer function for both mechanisms has been presented, and the region of highest receptivity of the film has been identified. In particular, while the liquid film behaves as a low pass filter against mechanisms B, a region of strong receptivity to mechanism A is found both upstream and downstream of the wiping point in the range of dimensionless frequencies 
$\hat {f}=0.03\text {--}0.05$. In galvanizing conditions, these correspond to wavelengths in the range 
$\lambda =25\text {--}35$ mm, in good agreement with industrial observations.
Finally, the comparison between the IBL and WIBL methods shows good qualitative agreement at the highest wiping numbers 
$\varPi _g$, highest rescaled Reynolds number 
$\delta$ and highest perturbation frequency 
$\hat {f}$. On the other hand, the TTBL method predicts a much larger thickness in the run-back flow. As mechanism A is triggered, this results in different wave shapes in the final coating. Nevertheless, these discrepancies do not alter the main results concerning the mechanisms on the undulation formation nor their range of large receptivity of the liquid film over the investigated wiping conditions. Although the interaction between the liquid film and impinging jet flow is characterized by coupling phenomena that do not fit in the simplified one-way coupling framework, the mechanisms revealed by this study certainly play an essential role in the stability of the jet wiping process.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2020.1075.
Acknowledgements
The authors gratefully acknowledge the essential contributions of M. Buszyk and D. Ninni. They took part in the development of the first IBL solver, implemented in the VKI software package BLEW (Boundary LayEr Wiping), during their short training program at the von Karman Institute. A.G. would like to thank CESGA (Centro de Supercomputación de Galicia) for computational resources.
Funding
The authors acknowledge financial support from ArcelorMittal. B.S. also thanks F.R.S.-FNRS for financial support.
Declaration of interest
The authors report no conflict of interest.
Appendix A. From Navier–Stokes equations to (1)
The long-wave formulation is based on the assumption that the streamwise reference scale 
$[x]$ is 
$[x]\gg [h]$ and, thus, a film parameter 
$\varepsilon =[h]/[x]\ll 1$ allows ordering of the importance of all the terms. The continuity equation, under the assumption of incompressibility, can be scaled as
\begin{equation} \left\{ \frac{U_p}{[x]}\right\} \partial_{\hat x} \hat{u}+\left\{ \frac{\varepsilon U_p }{[h]}\right\}\partial_{\hat y} \hat{v}=0, \end{equation}
and, thus, yields (2.1a) provided that 
$[v]=\varepsilon [u]=\varepsilon U_p$, having taken 
$[u]=U_p$. Taking 
$[t]=[x]/U_p$, the momentum equation in the streamwise coordinate can be scaled as
\begin{align} \left\{\frac{U_p^2}{[x]}\right\} (\partial_{\hat t}\hat{u}+ \hat {u} \partial_{\hat x} \hat{u} + \hat {v}\partial_{\hat y} \hat{u})=- \left\{ \frac{[p]}{\rho_l [x]}\right\} \partial_{\hat x} \hat{p}_l+\left\{\nu_l \frac{U_p}{[x]^2} \right\} \partial_{\hat x \hat x} \hat{u}+\left\{\nu_l \frac{U_p}{[h]^2} \right\} \partial_{\hat y\hat y} \hat{u} +g. \end{align}
Multiplying both sides by 
$[h]^2/\nu U_p$, taking the reference pressures 
$[p]$, thickness 
$[h]$ from table 1 and neglecting terms of 
${O}(\varepsilon ^2)$ yields (2.1b). The same procedure on the momentum equation along the cross-stream direction 
$y$ leads to (2.1c). The kinematic boundary conditions (2.2) can be obtained scaling by 
$\varepsilon U_p$, while the full force balance at the interface scales as
\begin{equation} (\hat p_l -\hat p_g) \{\rho_l g [x]\}- 2\mu_l \boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{E}_l\boldsymbol{\cdot} \boldsymbol{n}-2 \frac {[h] \sigma_l}{[x]^2} \hat{\kappa}=0 \quad\text{ in } y=h, \end{equation}
where the second term accounts for the viscous term in the normal direction and involves the normal unit vector 
$\boldsymbol {n}=(-\partial _{ x} {h},1)^T/\sqrt {1+(\partial _{x} h)^2}$ and the symmetric part of the rate of deformation tensor 
$\boldsymbol {E}_l=1/2\,(\boldsymbol {\nabla } \boldsymbol {v}+\boldsymbol {\nabla } \boldsymbol {v}^T)$, with 
$\boldsymbol {\nabla } \boldsymbol {v}=(u,v)$ the velocity field; the third term accounts for the surface tension, involving the mean curvature 
$\hat \kappa =-1/2\,\boldsymbol {\nabla }\boldsymbol {\cdot } \hat {\boldsymbol {n}}$. Expanding the viscous term in (A3) yields
\begin{equation} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{E}_l\boldsymbol{\cdot} \boldsymbol{n}=\frac{( \partial_y u + \partial_x v) \partial_x h- \partial_y v- \partial_x u (\partial_x h)^2 }{1+(\partial_x h)^2} \quad\mbox{ in } y=h. \end{equation} Scaling this expression and neglecting terms in 
${O}(\varepsilon ^2)$, the contribution of normal viscous stresses reads as
\begin{equation} 2\mu_l\,\frac{U_p}{[h]}{\hat {\boldsymbol{n}}\boldsymbol{\cdot}\hat{\boldsymbol{E}}_l\boldsymbol{\cdot}\hat{\boldsymbol{n}}}=2\mu_l \frac{U_p}{[h]} \varepsilon( \partial_{\hat y} \hat v-\partial_{\hat{x}}\hat{h} \partial_{\hat y}\hat u ) \quad\mbox{ in } \hat{y}=\hat{h}. \end{equation}
Introducing this result in (A3), observing that at 
${O}(\varepsilon )$ the dimensionless curvature becomes 
$\hat {\kappa }=1/2\,\partial _{\hat x\hat x} \hat h$ and dividing by the reference pressure 
$[p]=\rho \,g\,[x]$ gives
\begin{equation} (\hat p_l -\hat p_g)- \varepsilon^2{\hat{\boldsymbol{n}}\boldsymbol{\cdot}\hat{\boldsymbol{E}}_l \boldsymbol{\cdot}\hat{\boldsymbol{n}}}+\varepsilon^3\frac{l_c^2}{[h]^2} \partial_{\hat x \hat x}\hat{h}=0 \quad {\rm in}\ \hat y = \hat h, \end{equation}
where 
$l_c=\sqrt {\sigma /\rho \,g}$ is the capillary length. It is from this equation that the choice of the streamwise length scale 
$[x]$ – hence, the film parameter 
$\varepsilon$ – is taken. Following the scaling approach proposed by Shkadov (Reference Shkadov1971), this choice is made such that the surface tension contribution remains of leading order. Therefore,
\begin{equation} \varepsilon=\left(\frac{[h^2]}{l_c^2}\right)^{1/3}=\left(\frac{\mu U_p}{\sigma}\right)^{1/3}=Ca^{1/3}. \end{equation} By construction, then, the contribution of elongational viscosity is neglected and (A3) reduces to (2.3a). This approximation is valid for 
$Ca^{1/3}\ll 1$ and the weight of the viscous term becomes 
$Ca^{2/3}$. Finally, concerning the force balance in the tangential direction, the full equation reads as
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot}( 2\mu_l\boldsymbol{E}_l) \boldsymbol{\cdot}\boldsymbol{t} =\tau_g |_h,\text{ in } y=h, \end{equation}
where 
$\boldsymbol {t}=(1,\partial _x h)^T/\sqrt {1+(\partial _x h)^2}$ is the tangential unit vector. Expanding the matrix multiplication and dividing the result by 
$U_p/[h]$, the scaled form of this equation becomes
\begin{equation} (1-\varepsilon^2(\partial_{\hat{x}} \hat h)^2)(\varepsilon^2 \partial_{\hat{x}} \hat v +\partial_{\hat y} \hat u)+2 \varepsilon^2 \partial_{\hat x }\hat h (\partial_{\hat y}\hat {v}-\partial_{\hat x} \hat u )=\hat{\tau}_g \text{ in } \hat{y}=\hat{h}, \end{equation}
having introduced the reference shear stress 
$[\tau ]=\mu [u]/[h]$. To the leading order 
${O}(\varepsilon )$, this equation simplifies to (2.3b).
Appendix B. Details of the WIBL model derivation
After introducing the velocity profile (3.1), the fourth-order polynomial in 
$\bar {y}=y/h$ in (3.16) is obtained. Setting all the coefficients of the polynomial to zero the linear system in (3.17) is derived. The vector components on the right-hand side of (3.17), denoted as 
$\mathcal {G}=\{G_0,\ldots ,G_4\}$, are shown in table 5. Observe that all the variables in this and in the following tables are dimensionless and the ‘
$\hat{\,}$’ is dropped to ease the notation. Moreover, the derivation is carried out with an arbitrary scaling of the strip velocity, such that 
$u(y=0)=-\gamma$. This allows us to retrieve the model for the jet wiping if 
$\gamma =1$ and the classical models for falling liquid films if 
$\gamma =0$.
Table 5. Vector components 
$\mathcal {G}=\{G_0,\ldots ,G_4\}$ for the system in (3.17).
Setting all the coefficients 
$G_i$ equal to zero, the solution of the system in (3.17) yields the coefficients of the expansion 
$\mathcal {A}=\{a_0,\ldots a_4\}$ in table 6. These have a functional dependency on 
$a_0$ for all the terms weighted by the rescaled Reynolds number 
$\delta =\varepsilon Re$. Observe that in the falling film case, with 
$\gamma =0$ and 
$\partial _{x}p_g=\tau _g=0$, the equations in tables 5 and 6 recovers (6.43) and (6.44) in Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012.
Table 6. Full solution for the coefficients 
$\mathcal {A}=\{a_0,\ldots a_4\}$ in (3.1), obtained as vector 
$\mathcal {A}=\varGamma ^{-1} \mathcal {G}$ in (3.17).
The coefficients in table 6 are introduced in the contribution 
$\hat {q}_F$ of the flow rate (3.3), obtaining (B11) from table 7. From this, the first coefficient 
$a_0$ is isolated in (B12). Introducing this expression in the shear stress term (3.5) gives the results in (B13), having introduced 
$a_0= {3q}/{h}-{3}/{2}h\tau _g+3\gamma +{O}(\varepsilon )$ inside the parentheses and truncating all the terms in 
${O}(\varepsilon ^2)$. Equation (B13) (which, for 
$\gamma =1$, is (3.21)) can be finally introduced in (2.4b) to close the WIBL model.
Table 7. Full expression for 
$q$ using the coefficients 
$\mathcal {A}=\{a_0,\ldots a_4\}$ in table 6 in (3.3), together with the resulting expression of 
$a_0$ and the final result on the shear stress term 
${\rm \Delta} \tau$ for the WIBL model of the jet wiping process.
To retrieve the WIBL model for a falling liquid film (see Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012, (6.51)), it suffices to introduce 
$\gamma =0$ and 
$\partial _{x}p_g=\tau _g=0$ in (B12) and observe that in this case the advection term in (3.15) simplifies to 
$\mathcal {F}=6 q^2/(5h)$. Then, (2.4b) becomes (6.51) in Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012).
Appendix C. Implemented numerical schemes
The high-order fluxes 
$\boldsymbol {F}^{H}$ are derived from the two-step Lax–Wendroff scheme while the low-order fluxes 
$\boldsymbol {F}^{L}$ are taken from the two-step Lax–Friedrich scheme. The flux terms in the high-order scheme are
\begin{equation} \boldsymbol{F}^{H} = \begin{cases} \boldsymbol{F}^{+} = \boldsymbol{F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i+1}^{k})= \boldsymbol{F}({\boldsymbol V}_{i+{1}/{2}}^{k+{1}/{2}}), \\ \boldsymbol{F}^{-} = \boldsymbol{F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i-1}^{k})= \boldsymbol{F}({\boldsymbol V}_{i-{1}/{2}}^{k+{1}/{2}}). \end{cases} \end{equation}The flux terms in the lower-order scheme adds diffusive terms and reads as
\begin{equation} \boldsymbol{F}^{L} = \begin{cases} \displaystyle\boldsymbol{F}^{+} =\boldsymbol{F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i+1}^{k})=\frac 1 2 \boldsymbol{F}({\boldsymbol V}_{i+{1}/{2}}^{k+{1}/{2}})+\frac{{\rm \Delta} x}{2 {\rm \Delta} t} ({\boldsymbol V}_{i+1/2}^{k+1/2}-{\boldsymbol V}_{i}^{k}), \\ \displaystyle\boldsymbol{F}^{-} =\boldsymbol{F}({\boldsymbol V}_{i}^{k},{\boldsymbol V}_{i-1}^{k})=\frac 1 2 \boldsymbol{F}({\boldsymbol V}_{i-{1}/{2}}^{k+{1}/{2}})+\frac{{\rm \Delta} x}{2 {\rm \Delta} t} ({\boldsymbol V}_{i}^{k}-{\boldsymbol V}_{i-1/2}^{k+1/2}).\end{cases} \end{equation}Both schemes use the midpoint solutions in time:
\begin{gather} {\boldsymbol V}_{i-{1}/{2}}^{k+{1}/{2}} = \frac{1}{2}( {\boldsymbol V}_{i}^{k}+ {\boldsymbol V}_{i-1}^{k})+\frac{{\rm \Delta} t}{2{\rm \Delta} x}[ \boldsymbol{F}( {\boldsymbol V}_{i}^{k})-\boldsymbol{F}( {\boldsymbol V}_{i-1}^{k})] + \frac{1}{2}{\rm \Delta} t \boldsymbol{S}_{i+{1}/{2}}^{k}, \end{gather}
\begin{gather}{\boldsymbol V}_{i-{1}/{2}}^{k+{1}/{2}} = \frac{1}{2}( {\boldsymbol V}_{i}^{k}- {\boldsymbol V}_{i-1}^{k})+\frac{{\rm \Delta} t}{2{\rm \Delta} x}[ \mathbb{F}( {\boldsymbol V}_{i}^{k})-\boldsymbol{F}( {\boldsymbol V}_{i-1}^{k})] + \frac{1}{2}{\rm \Delta} t \boldsymbol{S}_{i-{1}/{2}}^{k}. \end{gather}
