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.

(3.6)\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

(3.7)\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

(3.8)\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

(3.9)\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:

(3.10)\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:

(3.11)\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

(3.12)\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

(3.13)\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

(3.14)\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

(3.15)\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}$:

(3.16)\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

(3.17)\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

(3.18)\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

(3.19)\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:

(3.20)\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

(3.21)\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.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

(4.1)\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}$,

(4.2)\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

(4.3)\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

(4.4)\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.

(4.5)\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:

(4.6)\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

(4.7)\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:

(4.8)\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

(4.9)\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):

(4.10)\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}

This term is closed using (4.7) with (4.3) and (4.2).

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

(6.1)\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

(6.2)\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

(6.3)\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

(6.4)\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:

(6.5)\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

(6.6)\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

(6.7)\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.