2.1 Governing equations

The steady Navier–Stokes equations describe the velocity $\boldsymbol{v}$ and pressure $p$ in the liquid phase:

The lengths, velocities and stresses are non-dimensionalized with the characteristic scales $Q/U$ , $U$ and $\unicode[STIX]{x1D707}_{adv}U^{2}/Q$ , respectively. The Reynolds, capillary and Bond numbers are given by $Re=\unicode[STIX]{x1D70C}Vw/\unicode[STIX]{x1D707}_{adv}$ , $Ca=\unicode[STIX]{x1D707}_{adv}U/\unicode[STIX]{x1D70E}$ and $Bo=(\unicode[STIX]{x1D70C}g/\unicode[STIX]{x1D70E})(wV/U)^{2}$ , respectively. The gravitational acceleration is given by $\boldsymbol{g}=-g\boldsymbol{e}_{y}$ , where $\boldsymbol{e}_{y}$ is the unit vector in the $y$ -direction.

Near the DCL, the air phase is long and slender. Lubrication theory is thus applied to the air phase, resulting in a 1D description of the air flow (as described in our previous work (Vandre Reference Vandre2013; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016)):

where $h$ is the interface height, $p$ is the air pressure, $u_{s}$ represents the horizontal interface velocity, $\unicode[STIX]{x1D706}=l_{slip}U/Q=l_{slip}/h_{inf}$ is a dimensionless slip length with $l_{slip}$ being the dimensional slip length, and $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D707}_{rec}/\unicode[STIX]{x1D707}_{adv}$ is the viscosity ratio. The velocity gradient in (2.3) is evaluated at the interface. Note that the quasi-parallel (QP) approximation ( $\text{d}x\approx \text{d}s$ ) is used in (2.2) and (2.3) to express pressure gradients with respect to the arclength coordinate $s$ along the interface (figure 2) (Jacqmin Reference Jacqmin2004; Sbragaglia, Sugiyama & Biferale Reference Sbragaglia, Sugiyama and Biferale2008).

The following set of dimensionless boundary conditions is applied along the upstream fluid interface:

Equations (2.4) and (2.5) represent the no-slip and no-penetration boundary conditions, and (2.6) and (2.7) are interfacial stress balances. At the fluid interface, $\boldsymbol{n}$ is the unit outward vector normal to the interface, and unit tangent vector $\boldsymbol{t}$ points in the direction of increasing distance $s$ along the interface (figure 2). Subscripts ‘ $adv$ ’ and ‘ $rec$ ’ indicate properties in the advancing and receding phases, respectively.

In the stress balances, $\unicode[STIX]{x1D735}_{\boldsymbol{s}}\unicode[STIX]{x1D70E}$ denotes the surface-tension gradient along the interface, $\unicode[STIX]{x1D705}$ represents the interface curvature, and $\unicode[STIX]{x1D64F}$ is the Newtonian stress tensor. The Marangoni and capillary numbers are defined as $M=(\unicode[STIX]{x1D70E}_{o}-\unicode[STIX]{x1D70E}_{m})/\unicode[STIX]{x1D707}_{adv}U$ and $Ca=\unicode[STIX]{x1D707}_{adv}U/\unicode[STIX]{x1D70E}_{m}$ . Here, $\unicode[STIX]{x1D70E}_{o}$ is the surface tension of the surfactant-free solution and $\unicode[STIX]{x1D70E}_{m}$ is the surface tension corresponding to the mean surfactant concentration $\unicode[STIX]{x1D6E4}_{m}$ at the interface. In the absence of surfactant, $\unicode[STIX]{x1D70E}_{o}=\unicode[STIX]{x1D70E}_{m}$ and the Marangoni stresses vanish (i.e. $M=0$ ).

At the downstream fluid interface, the no-slip and kinematic boundary conditions (2.4) and (2.5) are still valid. However, since the air flow near the downstream interface is not expected to have significant effects on wetting failure, the stress contribution from the air phase there is neglected. The stress balances at the downstream interface become:

where $P_{amb}$ is the ambient pressure. As noted earlier, we neglect the effect of surfactants on the downstream fluid interface. Therefore, surface tension there remains constant, resulting in (2.8) and (2.9).

To remove the stress singularity that would arise by applying the no-slip boundary condition at the DCL, a Navier slip boundary condition is applied along the moving substrate (Huh & Scriven Reference Huh and Scriven1971; Dussan Reference Dussan1976; Chan et al. Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013; Sibley, Nold & Kalliadasis Reference Sibley, Nold and Kalliadasis2015):

where $\unicode[STIX]{x1D706}$ is the dimensionless slip length defined earlier. Symbols $\boldsymbol{n}_{\boldsymbol{s}}$ , $\boldsymbol{t}_{\boldsymbol{s}}$ and $\boldsymbol{U}$ correspond to the substrate’s normal, tangent and velocity vectors, respectively. This slip boundary condition is applied along the entire bottom substrate and is able to recover the no-slip condition ( $\boldsymbol{v}\rightarrow \boldsymbol{U}$ ) for distances greater than $l_{slip}$ away from the DCL. Although the slip length will in general have different values in the advancing and receding phases, for simplicity we assume here that it has the same value in both phases. This assumption has worked well in our previous work on the parallel-plate system (Vandre et al. Reference Vandre, Carvalho and Kumar2012, Reference Vandre, Carvalho and Kumar2013, Reference Vandre, Carvalho and Kumar2014). Note that, in the vertical direction, the no-penetration condition is applied to the bottom boundary.

In the hydrodynamic model, the microscopic contact angle $\unicode[STIX]{x1D703}_{mic}$ serves as a boundary condition at the DCL to determine the free-surface shape (Blake Reference Blake2006; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). In general, $\unicode[STIX]{x1D703}_{mic}$ may depend on problem parameters such as the substrate speed, curtain height and feed flow rate (Kistler Reference Kistler and Berg1993; Shikhmurzaev Reference Shikhmurzaev1997; Sui et al. Reference Sui, Ding and Spelt2014; Blake et al. Reference Blake, Fernandez-Toledano, Doyen and De Coninck2015) (see § 1.1). Since we wish to isolate the influence of air stresses, we take the simplest view and assume $\unicode[STIX]{x1D703}_{mic}$ to be equal to the static contact angle $\unicode[STIX]{x1D703}_{s}$ (a fixed value) (Lowndes Reference Lowndes1980; Eggers Reference Eggers2005; Vandre et al. Reference Vandre, Carvalho and Kumar2013; Sui et al. Reference Sui, Ding and Spelt2014).

2.2 Surfactant transport

A steady convection–diffusion equation describes insoluble surfactant transport along the upstream fluid interface:

where $\unicode[STIX]{x1D735}_{\boldsymbol{s}}$ is the surface gradient operator and $\boldsymbol{v}_{\boldsymbol{s}}$ is the interface velocity. The surface Péclet number $Pe=Vw/D_{s}$ is defined using the surface diffusion coefficient $D_{s}$ . The dimensionless surfactant concentration $\unicode[STIX]{x1D6E4}$ is scaled with $K\unicode[STIX]{x1D6E4}_{m}$ , where $K$ represents a scaling factor ( $K\gg 1$ ) and $\unicode[STIX]{x1D6E4}_{m}$ is the mean surfactant concentration at the fluid interface. Because the flow along the interface carries surfactant towards the DCL, large concentration gradients arise there that are computationally prohibitive to resolve if we scale $\unicode[STIX]{x1D6E4}$ with $\unicode[STIX]{x1D6E4}_{m}$ only (i.e. $K=1$ ). As a result, we introduce a scaling factor $K$ and choose  $K=500$ .

The surfactant transport equation (2.11) is coupled with the tangential stress balance (2.6) through an equation of state (Schunk & Scriven Reference Schunk, Scriven, Kistler and Schweizer1997). For simplicity, we assume that the surfactant concentration is dilute enough to apply a linear equation of state (Ramé Reference Ramé2001; Campana et al. Reference Campana, Di Paolo and Saita2004; Kumar & Matar Reference Kumar and Matar2004; Campana et al. Reference Campana, Ubal, Giavedoni and Saita2011):

Here, surface tension $\unicode[STIX]{x1D70E}$ has been non-dimensionalized using the relation $\unicode[STIX]{x1D70E}^{\prime }=\unicode[STIX]{x1D70E}_{m}+(\unicode[STIX]{x1D70E}_{o}-\unicode[STIX]{x1D70E}_{m})\unicode[STIX]{x1D70E}$ , where $\unicode[STIX]{x1D70E}^{\prime }$ represents the dimensional surface tension (Craster & Matar Reference Craster and Matar2009).

We assume that surfactants do not deposit onto the solid surfaces. As a result, no-flux boundary conditions $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E4}=0$ are applied at the DCL and also the upstream tip of the slot die (figure 2). These boundary conditions require that a constraint of constant surfactant mass be applied:

where $S$ represents the total arclength of upstream fluid interface. Note that $1/K$ in (2.13) results from scaling concentration by $K\unicode[STIX]{x1D6E4}_{m}$ .

Surfactants can potentially influence the microscopic contact angle $\unicode[STIX]{x1D703}_{mic}$ and induce surface-tension gradients at the same time. To isolate the influence of Marangoni stresses on wetting failure, we assume that $\unicode[STIX]{x1D703}_{mic}$ is neutral at the DCL (i.e. $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ ) and remains unchanged by surfactants. In addition, when comparing curtain coating in the absence and presence of surfactants, the mean surface tensions $\unicode[STIX]{x1D70E}_{m}$ are kept the same.

2.3 Solution method

The Navier–Stokes equations (2.1a,b ), 1D description of the receding air flow (2.2), surfactant transport equation (2.11) and associated boundary and integral conditions (2.4)–(2.10), (2.12) and (2.13) are solved with the Galerkin finite element method. The liquid domain is partitioned into seven quadrangular regions and the air flow is only discretized in the 1D region (8) near the DCL (figure 3 a). Regions (8) and (9) represent the 1D domain on which surfactant transport along the upstream fluid interface is solved.

To solve the free-boundary problem, the physical domain (figure 3 a) is mapped to a computational domain whose coordinates are $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D709}$ (figure 3 b). The mapping function satisfies the elliptic mesh generation equations (Carvalho & Scriven Reference Carvalho and Scriven1997; Nam & Carvalho Reference Nam and Carvalho2009). Biquadratic basis functions are used for position and velocity variables, and linear discontinuous basis functions are used for pressure in the liquid phase. The air pressure and the surfactant concentration are also described by biquadratic basis functions.

Figure 3. Schematics of (a) the physical coating domain and (b) the computational domain composed of seven liquid regions (1)–(7), air region (8) and surfactant transport regions (8) and (9).

We found that mesh quality is best when the contact-line position $x_{DCL}$ (figure 3 a) is located directly underneath the upstream tip of the slot die ( $x_{DCL}=0$ ). When the DCL moves downstream ( $x_{DCL}>0$ ), elements near the DCL are stretched as the curtain is pulled (appendix A). In contrast, when $x_{DCL}<0$ , elements are stretched both near the DCL and near the upstream interface, where the curtain contacts a liquid heel that forms (appendix A). These mesh distortions make it very difficult to obtain discretization-independent steady-state solutions as the substrate speed $U$ increases (at a fixed feed flow velocity $V$ ). To mitigate this numerical issue, we fix $x_{DCL}$ (as a flow parameter) and calculate $V$ (as a dependent variable) as we trace solution paths to higher substrate speeds $U$ (using first-order continuation). As expected, $V$ must increase as $U$ increases to hold the DCL position fixed (appendix B).

Since the interface bends sharply near the DCL (figure 3 a), a larger number of elements near the DCL is required to resolve large gradients in the interface curvature. These gradients become larger as $Ca$ increases. As a consequence, we found that element sizes near the DCL need to be less than $10^{-3}l_{slip}$ to obtain mesh-independent solutions. Convergence is tested by systematically varying the number of elements to verify that solutions are mesh-independent (i.e. less than 2 % variation in $Ca^{crit}$ predicted from solution paths in § 3).

2.4 Model parameters

Because we are most interested in how flow behaviour changes as the substrate speed increases for fixed fluid properties, it is convenient to choose a representative set of dimensional parameters. Motivated by experimental data, we take $\unicode[STIX]{x1D707}_{adv}=25~\text{cP}$ , $\unicode[STIX]{x1D707}_{rec}=0.018~\text{cP}$ , $\unicode[STIX]{x1D70E}=70~\text{mN}~\text{m}^{-1}$ and $\unicode[STIX]{x1D70C}=1000~\text{kg}~\text{m}^{-3}$ . In addition, since we would like to isolate the effects of substrate speed and feed flow rate, we assume a neutral microscopic contact angle of $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ . We note that, even for this value of $\unicode[STIX]{x1D703}_{mic}$ , predictions of the hybrid model for the parallel-plate geometry are in excellent agreement with results from a model that uses a full 2D description of both the air and liquid phases (Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016). For this reason, we expect that a 1D description of the air flow is sufficient for accurately describing the behaviour of curtain coating.

To avoid the prohibitive computational costs that come with having too large a flow domain, we fix the curtain height at $h_{curt}=1~\text{cm}$ . Although this is smaller than typical curtain heights ( ${\sim}2{-}25~\text{cm}$ ), we expect that the results reported here will show the same qualitative behaviour as results for larger curtain heights.

Since smaller values of the slip length mean more elements near the DCL and a higher computational cost, a relatively large slip length $l_{slip}=10^{-5}~\text{m}$ is selected for numerical convenience. Other runs we have performed (not shown here) and our previous work with the parallel-plate geometry (Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016) indicate that an increase in the slip length only increases the $Ca^{crit}$ values without changing the qualitative nature of the solutions. We note that our choice of $l_{slip}$ leads to a situation where the slip length is larger than the thickness of the air film. However, this was also the case in our previous work on the parallel-plate system (Vandre et al. Reference Vandre, Carvalho and Kumar2013), yet the hydrodynamic model yielded predictions consistent with key experimental observations. This suggests that slip in the air phase, while necessary to include, may not play a dominant role in determining the main aspects of dynamic wetting failure for the flows considered in Vandre et al. (Reference Vandre, Carvalho and Kumar2013).

We consider values of $U$ (substrate speed) from ${\sim}0.1$ to $10~\text{m}~\text{s}^{-1}$ and $V$ (feed flow velocity) from ${\sim}0.1$ to $1~\text{m}~\text{s}^{-1}$ , and take the slot width $w$ to be 1 mm. With the parameter choices listed above, this yields values of $h_{inf}$ from ${\sim}0.1$ to 1 mm (appendix B). In addition, $Ca$ varies from ${\sim}0.1$ to 2.5, $Re$ varies from ${\sim}1$ to 40, $Bo$ varies from ${\sim}10^{-3}$ to $10^{-1}$ , and the dimensionless slip length $\unicode[STIX]{x1D706}$ varies from ${\sim}10^{-2}$ to $10^{-1}$ . Note that, since the liquid properties are fixed, $Ca$ , $Re$ , $Bo$ and $\unicode[STIX]{x1D706}$ change as the substrate speed and feed flow velocity are varied.

4.1 Flow fields with different $x_{DCL}$

Air/liquid displacement ( $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ ) in curtain coating is studied at various positions of $x_{DCL}$ , which is controlled by the feed flow velocity $V$ in experiments at fixed substrate speed $U$ and curtain height $h_{curt}$ (Blake et al. Reference Blake, Clarke and Ruschak1994; Miyamoto & Katagiri Reference Miyamoto, Katagiri, Kistler and Schweizer1997; Marston et al. Reference Marston, Simmons and Decent2007; Chang et al. Reference Chang, Shih, Liu and Tiu2012). Decreasing $V$ allows the substrate to drag the DCL downstream, whereas increasing $V$ tends to move the DCL upstream.

Figure 5 shows (a–c) pressure contours and (d–f) velocity-magnitude contours for three flow configurations (corresponding to three $x_{DCL}$ values) in curtain coating (Blake et al. Reference Blake, Clarke and Ruschak1994; Clarke Reference Clarke1995; Miyamoto & Katagiri Reference Miyamoto, Katagiri, Kistler and Schweizer1997; Marston et al. Reference Marston, Decent and Simmons2006; Chang et al. Reference Chang, Shih, Liu and Tiu2012): (a,d) bead pulling, $x_{DCL}>0$ ; (b,e) DCL right beneath the liquid curtain, $x_{DCL}=0$ ; and (c,f) heel formation, $x_{DCL}<0$ . Because it is easiest to compare different cases when the spatial scales are the same, all lengths have been rescaled by the slot width $w$ in this plot and subsequent ones. Note that the first two cases are at their respective critical capillary numbers $Ca^{crit}$ , whereas the last one is not. This is because, in the last case, mesh distortion is so strong that it prevents us from obtaining converged solutions before $Ca^{crit}$ can be clearly identified. Impinging liquid pressurizes the flow underneath the curtain, with the effect being most prominent when the DCL is right beneath the liquid curtain and in the heel-formation configuration. The effect is not as prominent in the bead-pulling configuration since the curtain is dragged by the substrate.

Figure 5. (a–c) Pressure contours and (d–f) velocity-magnitude contours for various $x_{DCL}$ values: $x_{DCL}=1$ (a,d), $x_{DCL}=0$ (b,e) and $x_{DCL}=-1$ (c,f). Values of other parameters are viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ . Here, $Ca=1.33$ ( $x_{DCL}=1$ ), 2.59 ( $x_{DCL}=0$ ) and 1.36 ( $x_{DCL}=-1$ ). Note that $Ca=Ca^{crit}$ except for the case $x_{DCL}=-1$ . In addition, $Re=11.99$  ( $x_{DCL}=1$ ), 30.37 ( $x_{DCL}=0$ ) and 34.93 ( $x_{DCL}=-1$ ), and $Bo=9.12\times 10^{-4}$ ( $x_{DCL}=1$ ),  $1.54\times 10^{-3}$ ( $x_{DCL}=0$ ) and $7.33\times 10^{-3}$ ( $x_{DCL}=-1$ ).

Prior experimental observations suggest that hydrodynamic assist is most influential (i.e. the substrate speed $U^{crit}$ is maximum) when the contact line is right beneath the liquid curtain (Blake et al. Reference Blake, Clarke and Ruschak1994, Reference Blake, Bracke and Shikhmurzaev1999; Yamamura et al. Reference Yamamura, Suematsu, Kajiwara and Adachi2000; Wilson et al. Reference Wilson, Summers, Shikhmurzaev, Clarke and Blake2006). Our results show that, when $x_{DCL}=0$ , $Ca^{crit}$ is the largest and its value is almost twice that when $x_{DCL}=1$ (figure 5 a,b). The degree of hydrodynamic assist is related to the distance between the DCL and the curtain impingement position (i.e. right beneath the upstream tip of the slot die) (figure 5 a–c). Hydrodynamic pressure due to the inertia of the liquid curtain has the strongest effect when the impingement position overlaps the DCL (i.e. $x_{DCL}=0$ ) because dynamic wetting and wetting failure take place at the DCL (i.e. air bubbles break off vertices of the sawtooth-shaped DCL) (Blake & Ruschak Reference Blake and Ruschak1979; Vandre et al. Reference Vandre, Carvalho and Kumar2014).

In contrast, in the bead-pulling configuration, the effect of hydrodynamic assist is weakened because the high-pressure region is not as prominent and the DCL is away from the impingement position (figure 5 a). In the heel-formation configuration, although hydrodynamic pressure is enhanced under the curtain, the increased distance between the impingement position and the DCL results in weak hydrodynamic assist (figure 5 c).

In all the flow configurations, liquid must accelerate near the DCL to match the substrate velocity (figure 5 d–f). Consequently, liquid pressure decreases sharply near the DCL, providing a pressure gradient that draws liquid towards the moving substrate (figure 5 a–c). The low-pressure zone paired with the high-pressure region under the curtain allow the $x_{DCL}=0$ case to achieve high coating speeds (figure 5 b). In the heel-formation configuration, however, the high pressure under the curtain is located far from the DCL and directs liquid upstream (figure 5 c). This results in a recirculation region (figure 5 f), similar to what is observed experimentally (Clarke & Stattersfield Reference Clarke and Stattersfield2006).

4.2 Mechanism of hydrodynamic assist

As noted in § 3, the mechanism of wetting failure is similar in curtain coating and the parallel-plate system. Here, we investigate the influence of hydrodynamic assist on the wetting failure mechanism. Figure 6(a) shows the capillary-stress gradient (filled symbols) and air-pressure gradient (open symbols) at the interface IP as a function of $Ca$ for $x_{DCL}=0$ (black circles) and $x_{DCL}=0.5$ (red diamonds). Notably, when $x_{DCL}=0$ the capillary-stress gradient is larger and has a stronger dependence on the IP location ( ${\sim}r_{f}^{-1.5}$ ) in comparison with the case when $x_{DCL}=0.5$ ( ${\sim}r_{f}^{-1}$ ). We point out that the $r_{f}^{-1}$ behaviour observed when $x_{DCL}=0.5$ is similar to that in the previously examined parallel-plate geometry (figure 4 b). Therefore, the stronger dependence of the capillary-stress gradient on $r_{f}$ when $x_{DCL}=0$ can be attributed to the presence of hydrodynamic assist.

Figure 6. (a) Magnitude of capillary-stress (filled symbols) and air-pressure (open symbols) gradients at the IP for $x_{DCL}=0$ (black circles) and $x_{DCL}=0.5$ (red diamonds). (b) IP radial distance $r_{f}$ (filled symbols) and height $h_{f}$ (open symbols) as a function of $Ca$ for $x_{DCL}=0$ (black circles) and $x_{DCL}=0.5$ (red diamonds). Inset: the ratio between the two length scales for $x_{DCL}=0$ (black circles) and $x_{DCL}=0.5$ (red diamonds). Values of other parameters are viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ .

Figure 6(b) shows the change of the interface IP height $h_{f}$ (filled symbols) and radial distance $r_{f}$ (open symbols) with $Ca$ (see figure 2 for definitions of $h_{f}$ and $r_{f}$ ). For both $x_{DCL}=0$ (black circles) and $x_{DCL}=0.5$ (red diamonds), the IP migrates towards the DCL (i.e. $h_{f}$ and $r_{f}$ decreasing) as $Ca$ increases because the interface curvature needs to increase to balance the increasing viscous stresses from the air. Notably, the IP for $x_{DCL}=0$ is closer to the DCL than that for $x_{DCL}=0.5$ as $Ca$ approaches the onset of wetting failure, meaning that the interface curvature for $x_{DCL}=0$ is larger than that for $x_{DCL}=0.5$ . As a result, the capillary-stress gradient for $x_{DCL}=0$ is larger than that for $x_{DCL}=0.5$ , which is consistent with what is shown in figure 6(a).

The inset of figure 6(b) shows the ratio between $h_{f}$ and $r_{f}$ for $x_{DCL}=0$ and $x_{DCL}=0.5$ . This ratio is related to the angle of the wedge formed in the receding phase (figure 2); a larger ratio corresponds to a smaller macroscopic contact angle $\unicode[STIX]{x1D703}_{M}$ . The inset of figure 6(b) shows that this ratio is larger when $x_{DCL}=0$ . This is a consequence of the strong hydrodynamic pressure near the DCL when $x_{DCL}=0$ , which suppresses the penetration of the air film. Therefore, the model results suggest that the macroscopic contact angle $\unicode[STIX]{x1D703}_{M}$ in curtain coating is reduced by hydrodynamic assist, which is consistent with experimental observations that hydrodynamic assist leads to a reduction in the dynamic contact angle $\unicode[STIX]{x1D703}_{D}$ (Blake et al. Reference Blake, Bracke and Shikhmurzaev1999, Reference Blake, Dobson and Ruschak2004; Clarke & Stattersfield Reference Clarke and Stattersfield2006; Decent Reference Decent2008). In § 5, we will discuss the influence of air stresses on the contact-angle behaviour.

Together, figures 5 and 6 suggest that hydrodynamic assist can be attributed to the strong pressure generated by the inertia of the impinging curtain. This pressure creates large gradients in the interface curvature near the DCL, which in turn enhances the capillary-stress gradient there and lowers the dynamic contact angle. These larger capillary-stress gradients more effectively pump air away from the DCL and delay the onset of wetting failure.

6.1 Solution paths

The solution paths for different Marangoni numbers when $x_{DCL}=0$ are shown in figure 9. We choose to fix $x_{DCL}=0$ because this is where hydrodynamic assist is maximum, in both the presence and the absence of surfactants. The solution path for $M=0$ (circles) corresponds to the case where surfactants are absent (figure 7). Notably, the presence of Marangoni stresses ( $M>0$ ) decreases $Ca^{crit}$ and promotes the onset of wetting failure, suggesting a possible mechanism for the experimental observations of Marston et al. (Reference Marston, Hawkins, Decent and Simmons2009). Figure 9 also demonstrates that increasing the strength of Marangoni stresses (i.e. higher $M$ or surfactant concentration) will further decrease $Ca^{crit}$ , which was also experimentally observed by Marston et al. (Reference Marston, Hawkins, Decent and Simmons2009).

Figure 9. Solution paths for various $M$ when $x_{DCL}=0$ : $M=0$ (black circles), $M=0.15$ (red diamonds), $M=0.3$ (blue squares) and $M=0.75$ (magenta triangles). Inset: solution paths for $M=0$ (black circles) and $M=0.15$ (red diamonds) when $x_{DCL}=0.5$ . Values of other parameters are viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ . For $x_{DCL}=0$ , $Re=30.37$ and $Bo=1.54\times 10^{-3}$ at $Ca=Ca^{crit}$ when $M=0$ . For $x_{DCL}=0.5$ , $Re=15.07$ and $Bo=7.94\times 10^{-4}$ at $Ca=Ca^{crit}$ when $M=0$ .

The inset of figure 9 shows the solution paths for $M=0$ and $M=0.15$ at $x_{DCL}=0.5$ . As noted in § 5, the feed flow rate at $Ca^{crit}$ decreases as the DCL moves downstream. Therefore, the feed flow rates for solution paths in the inset of figure 9 are smaller than those leading to $x_{DCL}=0$ in figure 9. Notably, in going from $x_{DCL}=0$ to $x_{DCL}=0.5$ , $Ca^{crit}$ decreases 30 % for the surfactant-free case (i.e. $M=0$ ). In contrast, the reduction in $Ca^{crit}$ for $M=0.15$ is only 6.3 %, suggesting that the influence of the feed flow rate on the critical substrate speeds (i.e. the degree of hydrodynamic assist) is weakened by Marangoni stresses. This decreased degree of hydrodynamic assist is similar to what is observed by Marston et al. (Reference Marston, Hawkins, Decent and Simmons2009), again indicating that Marangoni stresses may play an important role in curtain coating.

6.2 Mechanism of influence of surfactants

Our previous work on the parallel-plate system (figure 1 b) suggests that Marangoni stresses promote the onset of wetting failure by thinning the air film between the fluid interface and the substrate (Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016). However, since curtain coating has a significantly different geometry, it is not obvious that it shares the same physical mechanism. To address this issue, we perform a stress-gradient analysis.

Figure 10. (a) Magnitude of stress gradients at the interface IP for $M=0$ (black circles) and $M=0.15$ (red diamonds). The capillary-stress gradients (filled symbols) match the air-pressure gradients (open symbols) at $Ca^{crit}$ values denoted by dotted lines. (b) IP radial distance $r_{f}$ (filled symbols) and height $h_{f}$ (open symbols) as a function of $Ca$ for $M=0$ (black circles) and $M=0.15$ (red diamonds). Inset: the ratio between the two length scales for $M=0$ (black circles) and $M=0.15$ (red diamonds). Values of other parameters are $x_{DCL}=0$ , viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ .

Capillary-stress gradients (filled symbols) and air-pressure gradients (open symbols) at the interface IP are plotted as a function of $Ca$ in figure 10(a) when Marangoni stresses are absent ( $M=0$ , black circles) and present ( $M=0.15$ , red diamonds). Similar to what we found in our previous work on the parallel-plate geometry (Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016), the presence of Marangoni stresses tends to increase the air-pressure gradient (red open diamonds) such that it matches the capillary-stress gradient (red filled diamonds) at a smaller $Ca^{crit}$ compared to when  $M=0$ .

Flow along the interface carries surfactant towards the DCL, creating large concentration gradients there (figure 11). The surfactant concentration is largest near the DCL, leading to a Marangoni stress directed away from the DCL. Stronger shear stresses in the air phase are required to balance Marangoni stresses along the fluid interface (see (2.6)), and this requires air velocity changes over a smaller distance ( $\text{d}u/\text{d}y$ ) near the DCL. Therefore, the air film near the DCL is thinned (i.e. $h_{f}$ decreases), increasing the magnitude of the air-pressure gradient, which scales as $h_{f}^{-2}$ (Vandre et al. Reference Vandre, Carvalho and Kumar2013; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016).

Figure 11. Surfactant concentration profile for $M=0.15$ and $Ca=0.93$ . The DCL is at $y=0$ . Values of other parameters are $x_{DCL}=0$ , viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ .

Figure 10(b) demonstrates the migration of the IP, where the characteristic length scales ( $r_{f}$ , filled symbols; $h_{f}$ , open symbols) are plotted as a function of $Ca$ . Notably, both $r_{f}$ and $h_{f}$ are smaller when $M=0.15$ (diamonds) compared to when $M=0$ (circles), indicating that the size of the air film shrinks in the presence of Marangoni stresses. Although smaller $r_{f}$ values indicate a larger interface curvature and a stronger capillary-stress gradient, the decrease of $h_{f}$ (with $Ca$ ) when $M=0.15$ is so drastic that the air-pressure gradient increases rapidly and becomes equal to the capillary-stress gradient at a smaller $Ca$ than that when  $M=0$ .

The ratio between the two length scales ( $h_{f}/r_{f}$ ) for $M=0$ (circles) and $M=0.15$ (diamonds) is shown in the inset of figure 10(b). Owing to the large reduction in $h_{f}$ , the ratio for $M=0.15$ is much smaller than that for $M=0$ . A smaller ratio corresponds to a larger macroscopic contact angle $\unicode[STIX]{x1D703}_{M}$ at the IP (figure 2), which is consistent with the larger values of $\unicode[STIX]{x1D703}_{M}$ in the solution paths when $M>0$ (figure 9).

As discussed in § 4, when Marangoni stresses are absent ( $M=0$ ), the strong dependence of the capillary-stress gradient on $r_{f}$ when $x_{DCL}=0$ can be attributed to the large pressure changes near the DCL ( ${\sim}r_{f}^{-1.5}$ in figure 10 a). When Marangoni stresses are present ( $M=0.15$ ), figure 10(a) shows that the dependence of capillary-stress gradients on $r_{f}$ is weakened ( ${\sim}r_{f}^{-1.1}$ ). To understand the reason for this, flow fields at the onset of wetting failure (i.e. at the respective $Ca^{crit}$ ) for $M=0$ and $M=0.15$ are presented in figure 12.

Figure 12 shows (a–c) pressure contours and (d–f) velocity-magnitude contours for $M=0$ and $M=0.15$ . Panels (a,b) and (d,e) are at the onset of wetting failure, and we note that the value of $Ca^{crit}$ when $M=0$ is larger than that when $M=0.15$ . Thus, panels (c,f) are for $M=0$ at a value of $Ca$ comparable to $Ca^{crit}$ for $M=0.15$ . In all cases, $x_{DCL}=0$ .

Figure 12. (a–c) Pressure contours and (d–f) velocity-magnitude contours for $M=0$ (a,c,d,f) and $M=0.15$ (b,e) at the onset of wetting failure. Values of other parameters are $x_{DCL}=0$ , viscosity ratio $\unicode[STIX]{x1D712}=7.2\times 10^{-4}$ and $\unicode[STIX]{x1D703}_{mic}=90^{\circ }$ . (a,d)  $Ca=2.59$ , $Re=30.37$ and $Bo=1.54\times 10^{-3}$ . (b,e)  $Ca=0.93$ , $Re=21.41$ and $Bo=5.89\times 10^{-3}$ . (c,f)  $Ca=0.96$ , $Re=20.94$ and $Bo=5.3\times 10^{-3}$ .

When $M=0$ (figure 12 a,d), the curtain remains nearly perpendicular to the substrate and the impinging liquid pressurizes the flow underneath the curtain. In contrast, a liquid bulge is formed at the upstream fluid interface near the DCL when $M=0.15$ (figure 12 b,e). This liquid bulge results from the presence of Marangoni stresses, which counteract the impinging liquid and slow down the liquid velocity near the fluid interface and the DCL (figure 12 e). We note that the bulge is not as prominent when $M=0$ and $Ca$ has a comparable value to $Ca^{crit}$ for $M=0.15$ (figure 12 c,f). The low-speed region when $M=0.15$ (figure 12 e) indicates that inertial forces near the DCL are weakened, suggesting a weaker influence of the impinging liquid on the DCL compared to when $M=0$ . This weakened influence may explain the weaker dependence of the capillary-stress gradient on $r_{f}$ when $M=0.15$ (figure 10 a) and also the decrease in the degree of hydrodynamic assist observed by Marston et al. (Reference Marston, Hawkins, Decent and Simmons2009).

Finally, we comment on several other phenomena that may be present in the experiments but are not accounted for in our calculations. First, the surfactants used in the experiments of Marston et al. (Reference Marston, Hawkins, Decent and Simmons2009) were soluble, whereas our model assumes that the surfactants are insoluble. We expect that surfactant solubility will weaken Marangoni effects since surfactants could desorb from the fluid interface into the bulk. Second, our model does not account for the possibility that surfactants may adsorb to the substrate. Third, the introduction of surfactants will lower the mean surface tension and the microscopic contact angle, which were taken as constant in the above calculations to isolate the influence of Marangoni stresses. Lowering $\unicode[STIX]{x1D703}_{mic}$ would be expected to raise the value of $Ca^{crit}$ since the substrate would be more wettable. As this is the opposite of what is observed in the experiments, it is not likely to be the cause of the decrease in $Ca^{crit}$ . Accounting for these phenomena as well as larger curtain heights will be important to more fully understand the experimental observation that surfactants significantly promote the onset of wetting failure.