*Article* **Coating Flow Near Channel Exit. A Theoretical Perspective**

## **Roger E. Khayat \* and Mohammad Tanvir Hossain**

Department of Mechanical and Materials Engineering, University of Western Ontario, London, ON N6A 5B9, Canada; mhossa69@uwo.ca

**\*** Correspondence: rkhayat@uwo.ca; Tel.: +1-519-661-2111

Received: 13 September 2020; Accepted: 9 October 2020; Published: 15 October 2020

**Abstract:** The planar flow of a steady moving-wall free-surface jet is examined theoretically for moderate inertia and surface tension. The method of matched asymptotic expansion and singular perturbation is used to explore the rich dynamics near the stress singularity. A thin-film approach is also proposed to capture the flow further downstream where the flow becomes of the boundary-layer type. We exploit the similarity character of the flow to circumvent the presence of the singularity. The study is of close relevance to slot and blade coating. The jet is found to always contract near the channel exit, but presents a mild expansion further downstream for a thick coating film. We predict that separation occurs upstream of the exit for slot coating, essentially for any coating thickness near the moving substrate, and for a thin film near the die. For capillary number of order one, the jet profile is not affected by surface tension but the normal stress along the free surface exhibits a maximum that strengthens with surface tension. In contrast to existing numerical findings, we predict the existence of upstream influence as indicated by the nonlinear pressure dependence on upstream distance and the pressure undershoot (overshoot) in blade (slot) coating at the exit.

**Keywords:** coating flow; free surface; boundary layer; stress singularity; matched asymptotic expansions

## **1. Introduction**

We examine the free-surface flow of a planar moving-wall jet at moderate Reynolds and capillary numbers near and far from the channel exit as encountered in coating flow. With the advent of high-speed coating and the use of low-viscosity liquids, inertia is becoming increasingly important but has been traditionally neglected in the modelling of coating flow. General aspects on the classification and analyses of coating flows can be found in the reviews by Ruschak [1] and Weinstein and Ruschak [2]. The current work focuses on slot and blade coating flows which, in addition to the substrate movement, typically involve an adverse or a favorable pressure gradient.

In slot coating, the liquid is forced into the slot die, and distributes through the narrow slot before it emerges onto the moving substrate. A low-pressure area or vacuum is imposed upstream of the die to facilitate a faster and more stable coating process. Consequently, part of the flow is driven downstream by the moving substrate and part of it circulates upstream in the low-pressure area, which causes a streamwise adverse pressure gradient to act inside the channel formed between the downstream die and the moving substrate. In blade coating, the flow lies between a fixed blade of a prescribed shape and a substrate moving parallel to itself. The coating liquid is dragged inside the channel by the moving substrate, which causes a hydrodynamic pressure rise at the upstream of the blade. The pressure rise causes the blade to reject most of the liquid, and only a fraction passes into the narrow channel. Since the drag flow can only carry half of the coating liquid, a streamwise favorable pressure gradient is generated to carry the rest [3]. In this case, the pressure gradient forces the coating liquid in the same direction as the movement of the substrate inside the channel. In both slot and blade

coating, the moving substrate drags the flow out of the channel in the form of a free-surface wall jet, and a thin layer of liquid film is obtained. In the present work, only the pressure gradient immediately upstream of the blade or die exit is accounted for, which is found in reality to be very close to constant. Consequently, the flow inside the channel will be assumed to be a superposition of Couette (velocity driven) and Poiseuille (pressure driven) flows.

Aside from its industrial importance, coating flow is fundamentally important as it spans various flow regimes. In an effort to clearly identify these regimes, de Ryck and Quere [4] carried out different experiments on fibre coating and used dimensional arguments, which illustrate the situation of coating flow in general. The visco-capillary range corresponds to slow coating flow and negligible Weber number, obeying the well-known Landau equation [5]. In this case, the film thickness results from a balance between capillarity and viscous forces. For coating flow of liquids of low viscosity (such as water) at higher velocity (about 1 m/s), the measurements of de Ryck and Quere were shown to deviate considerably from the Landau law, even if the capillary number remained negligible (Ca < 0.05). This is the visco-inertial regime, corresponding to a Weber number of order 1. At larger velocities (We > 1), the coating flow does not seem to depend any longer on the surface tension. This observation corroborates well findings for coating flow in general, as the thickness measurements reported by Lee et al. [6] and Becker and Wang [7] indicate for slot coating. De Ryck and Quere [4] also point out that the boundary-layer regime should be relevant to most industrial fibre-coating processes. Similarly, the more recent thickness measurements of Chang et al. [8] for low-viscosity slot-coating films show that viscous and surface tension effects become important only below a critical Reynolds number. At low speed, the minimum wet thickness increases with increasing capillary number but becomes independent of the capillary number for Ca > 0.3. Above the critical Reynolds number, fluid inertia becomes dominant. In this region, the minimum wet thickness decreases as the Reynolds number increases. Chang et al. [8] also established that inertial forces cannot be neglected in slot coating for Re larger than 10, which roughly corresponds to a substrate speed of 1 m/s. The earlier measurements of Carvalho and Kheshgi [9] show that the visco-capillary model remains valid at capillary numbers below unity. Their measurements also indicate that the model becomes less and less useful as the capillary number and the Reynolds number rise beyond unity. According to Romero et al. [10], high capillary number (Ca ≈ 5) and Reynolds number (Re ≈ 3) occur in higher-speed coating operations. It is this inertial regime that is of primary interest in the present study.

Although inertia has the desirable effect of lowering the flow limit and stabilizing the coating process [9], existing studies have predominantly focused on the role of surface tension. Lee et al. [6] and Chin et al. [11] showed that reducing the viscosity lowers the minimum thickness and increases the maximum coating speed. They found that a higher viscosity of the coating solution tends to be destabilizing, and a lower coating speed under high viscosity induces coating defects such as air entrainment and ribbing. Lee et al. [6] measured the minimum wet thickness for extrusion slot coating. In their experiment, the range of Reynolds number was low, from 0.2 to 24. They noticed that there exists a critical capillary number beyond which the minimum wet thickness becomes constant regardless of the value of Ca. Following the work of Lee et al. [6], and using a highly viscous fluid in the slot coater, Yu et al. [12] found that a minimum wet thickness can be reached if a low viscous fluid is used as a bottom carrier layer. Chang et al. [8] visualized the slot coating process for a low viscosity fluid, and observed that beyond a critical Reynolds number (Re = 20), both viscous and surface tension effects become negligible, with inertia dominating the flow.

Although extensive theoretical work has been devoted to the analysis of coating flow, the focus has been primarily on inertialess flow of a non-Newtonian fluid [13]. To a much lesser extent, numerical studies accounting for inertia were also carried out. An overview of high-speed blade coating was given by Aidun [14], with direct relevance to paper coating. One of the earlier numerical work was done by Saito and Scriven [15], who carried out a finite-element analysis coupled to an iterative scheme to examine the capillary number effect on the curved meniscus close to the static contact line in slot coating. They showed that the downstream meniscus no longer remains attached to the slot die at higher capillary number and lower flow rate. They found that the rate of convergence is highly dependent on the capillary number, and became exceedingly slow for Ca > 10. Carvalho and Kheshgi [9] and Jang and Song [16] examined the low-flow limit for slot coating in the inertial-capillary regime. Iliopoulos and Scriven [17] examined numerically the influence of particle suspensions on the blade and substrate. They found that the blade coating thickness increased with the Reynolds number, for any capillary number, but the coating thickness tapered off at the higher capillary numbers. In their study, the approximate ranges were, 0 < Re < 55 and 10 < Ca < 30. Lin et al. [18] carried out a comparison between numerical predictions and experiment, assessing operating windows for slot coating for their Reynolds number ranging from 0 to 100.

Perhaps the scarcity of theoretical studies on inertial coating flow is due to the stronger singularity at the contact point (line) at higher Reynolds number. Bajaj et al. [19] reported that for surface-tension dominated flow, a circulation appears upstream of the exit. This recirculation zone was found to gradually shrink with increasing Reynolds number. Consequently, a geometric singularity emerges in slot coating flow, resulting from the discontinuous slope at the contact line, which is increasingly exposed in relatively stronger flows. The singularity is also physical as a result of the change in the boundary condition from adherence at the die wall to the shear-free condition on the free surface (see also below). Inertia can be used to counteract the receding action of the downstream meniscus and the flow reversal, and delay the onset of the flow limit [9]. However, this gain will be offset by the strengthening of the singularity, which is no longer attenuated by a dominant surface tension. In this case, the circulation just upstream of the contact line weakens and eventually disappears altogether with increasing inertia, causing a more abrupt drop in the shear stress and a stronger singularity. We recall that the flow limit corresponds to a minimum coating thickness for a substrate speed or a maximum substrate speed for a given coating thickness. Inertia helps reduce the flow limit when, for instance, lower viscosity liquids are used, which tend to spread more easily [6]. The measurements and theoretical predictions of Carvalho and Keshghi [9] indicate that thinner films can be obtained at faster web speeds. Finally, the imposition of appropriate boundary conditions for free-surface inlet and outlet flows remains an open issue, particularly for small or vanishing surface tension [20].

Coating flow is modelled as a free-surface wall jet emerging from a channel bounded by a semi-infinite stationary (die or blade) wall and an infinite moving substrate. We assume inertia to remain relatively important, allowing the asymptotic development of the flow field in terms of some inverse power of the Reynolds number. The jet is assumed to be subject to a constant favorable or adverse pressure gradient far upstream where the flow acquires a Couette–Poiseuille (CP) character, typically as encountered in blade and slot coating. The assumption of fully developed upstream flow and the channel flow geometry are generally reasonable and commonly adopted in numerical simulation [15,19,21,22]. In roll coating, for instance, it is established that when the radius of the roller is large compared to the capillary length <sup>q</sup> σ ρg , then the effect of the curvature of the cylinder is very weak and results in only small perturbations to the uniform film thickness on the roller [23]. Here σ is the surface tension coefficient, ρ is the fluid density and g is the gravitational acceleration.

The stress singularity constitutes the major difficulty in a computational approach. The incorporation of the singularity point and its immediate vicinity is unavoidable in this case since the entire flow domain must be considered (discretized). The neighborhood region around the singularity, which is crucial to the rest of the flow domain, is difficult to handle numerically. In a numerical approach, the singularity is typically smoothed out or smeared over. The present asymptotic approach represents a viable alternative, at least for the flow in the vicinity of the singularity, which is avoided altogether as a result of the flow similarity in the free-surface layer and the boundary layer near the wall. In a numerical approach, the mesh is typically refined near the singularity, thus capturing more closely the singular behavior; the numerical difficulty resides in handling the resulting stronger flow gradients. Mitsoulis [24] conducted a finite-element analysis of blade coating. From his figures, we can see that both the shear stress as well as the pressure and the normal stress become singular at the channel exit. Finally, although the flow is primarily dominated by the stress singularity close

to the exit, this eventually influences the accuracy of the numerical predictions of the flow further downstream and upstream.

Saito and Scriven [15] used an iterative numerical scheme to determine the free surface profile. They reported on earlier studies [25] where convergence difficulties were encountered. This was particularly the case for slot coating, with severe bending of the free surface. Convergence difficulties were also encountered when inertia grows large at a fixed flow rate and dominates surface tension. Saito and Scriven [15] were able to overcome the convergence difficulties by combining simple coordinate parametrizations of parts of the free surface, and computing accurately the derivatives of finite-element residuals with respect to free-surface locations along the spines. Alternatively, mapping techniques have also been used more recently, where the governing equations and boundary conditions are transformed to an equivalent set defined in a known reference frame [9,21,22]. The volume-of-fluid approach has also been employed [16].

The present asymptotic approach circumvents the singularity altogether and does not require an iterative scheme to determine the shape of the meniscus. Perhaps more importantly, the present formulation provides a deeper insight on the flow structure and the rich dynamics near the singularity (see Section 7 for further discussion on this point and the relation with the triple-deck approach). We observe that the use of the singular perturbation technique also provides a procedure by which higher-order terms could be determined to the desired accuracy [26]. The formulation is essentially analytical, providing the steady flow solution in a manner that is completely amenable to the implementation of a linear or nonlinear stability analysis. Generally, asymptotic analyses have been successfully adopted for flows in the visco-capillary range [27,28]. Timoshin [29] developed a high-Reynolds-number asymptotic theory and examined the stability of boundary-layer flow over a coated surface. More recently, Tsang et al. [30] implemented a high-Reynolds-number asymptotic approach to study the interaction of a boundary layer on a solid plate and the free surface above. Studies of closer relevance to the present flow in the visco-inertial range were also conducted, but to a much lesser extent. Tillett [31] analyzed the laminar symmetric free-jet flow near the channel exit using the method of matched asymptotic expansions. Miyake et al. [32] carried out a similar analysis on a vertical jet of inviscid fluid, taking into account gravity effect. Philippe and Dumargue [33] applied an analysis similar to Tillett's for viscous axisymmetric vertical jets, emphasizing the interplay between the effects of gravity and inertia on the free surface shape and the velocity profile. Their approach, which is similar to the one in the present study, was validated against experiment. A local similarity transformation was carried out by Wilson [34] for the axisymmetric viscous-gravity jet emerging from a tube, but, unlike the present formulation, Wilson neglected the upstream influence. Khayat and co-workers examined the flow of a jet emerging from a channel of Newtonian [35,36] and non-Newtonian [37,38] liquids. We refer the reader to the book by Sobey [39] on interactive boundary layer and triple-deck theory for a perspective on asymptotic analyses, their applications and historic development.

Asymptotic analysis has also been applied to study coating flow, but the majority of the studies focused on dominant surface tension. Ruschak [26] performed a theoretical analysis based on the thin-film theory of Landau and Levich [5] to determine the effects of different parameters on the flow limit of extrusion slot coating for a Newtonian fluid. Ruschak considered a very small capillary number by setting the coating speed close to zero, and determined the film thickness in the slow flow (negligible inertia) limit by carrying out a singular perturbation method. Ruschak observed that film thickness becomes thinner for smaller capillary number, and noticed that gravity does not affect the flow much in the limit considered in the study. Higgins and Scriven [40] extended Ruschak's analysis by incorporating the viscous effect in the coating bead (inside the channel) with variable meniscus location. They concluded that as the coating speed increases, the dynamic contact angle increases, resulting in an altered coating behavior. Christodoulou and Scriven [41] provided an asymptotic analysis based on a thin-film approach in slide coating. Carvalho and Kheshgi [9] assessed the low-flow limit for slot coating, which is defined as the maximum substrate speed possible without any coating defects at a fixed minimum thickness or vice versa. By modifying the visco-capillary model, Carvalho and Kheshgi [9] developed a 2D numerical tool to study the effect of higher capillary numbers. They observed that at higher capillary number (Ca > 0.3), inertia started to dominate the flow, delaying the onset of the low-flow limit. Of close relevance to the present formulation is the matched asymptotic approach used by Rushack and Scriven [42], focusing on low flow rate and high surface-tension limits for an inertialess jet.

Other coating configurations were also analyzed. Savage [43] and, later, Gaskell et al. [44] analyzed the meniscus roll-coating flow of an inertialess fluid. Kelmanson [45] examined the effect of inertia in the small- and large- surface-tension limits for a coating film flowing around a rotating cylinder. Blythe and Simpkins [46] developed an asymptotic analysis using the Landau–Levich equation to examine the inertialess falling coating flow on a fibre. Later, Jang et al. [16] proposed a model for slot coating by modifying the visco-capillary model to accommodate the high capillary and to some extent inertia effects. Unlike the visco-capillary model, they accounted for the pressure variation under the slot die, which enabled them to predict the coating thickness at high capillary and Reynolds numbers. They observed that when Re > 10 the film thickness decreased with increasing inertia, and visco-capillary model did not show good agreement with the simulation for Ca > 0.2.

Finally, the present work is reminiscent of the early treatment of Goldstein [47] of the flow near the trailing edge of a semi-infinite plate, where the flow downstream of the singular edge was examined. Later, Stewrtson [48] and Messiter [49] developed a more comprehensive treatment of the flow very near the singularity, laying down the foundation for the triple-deck theory [39]. They established a rational theoretical approach for the investigation of separation and other non-linear features emerging in external flow situations. The present wall-jet flow differs fundamentally from the external separation flow due to the absence of an external inviscid flow. The present development is somewhat of the same level as the asymptotic treatment of Goldstein. The more complete analysis of the full Navier–Stokes equations and the implementation of a triple-deck approach may be envisageable in a future extension.

#### **2. Problem Formulation and the Physical Domain**

We consider the planar flow of an incompressible Newtonian fluid of density ρ, viscosity µ and surface tension σ, flowing between the semi-infinite stationary blade or slot die and the translating substrate, separated by a distance D. The flow is induced by the simultaneous forward translation velocity V of the substrate and an applied adverse or favorable pressure gradient g far upstream, where fully-developed Couette–Poiseuille (CP) flow is assumed to prevail. The assumption of fully-developed flow far upstream is commonly adopted in numerical studies on slot coating [10,15,21]. We also assume the CP flow to hold for blade coating as the pressure gradient is sensibly constant over a sufficiently long upstream distance [17]. Far downstream, the flow becomes uniform with speed V and film thickness T.

Non-dimensional variables are introduced by measuring the coordinates (x, z), the velocity components (u, w), the stream function ψ and the pressure p in units of D, V, VD and ρ*V* 2 , respectively. There result three dimensionless parameters appearing in the problem, namely the Reynolds number Re, the capillary number Ca, and the thickness-to-gap ratio Q, which are expressed as

$$\text{Re} = \frac{\rho \text{VD}}{\mu}, \qquad \text{Ca} = \frac{\mu \text{V}}{\sigma}, \qquad \text{Q} = \frac{\text{T}}{\text{D}}.\tag{1}$$

A parameter related to Q, which is conveniently introduced here:

$$\mathbf{G} \equiv \mathbf{3} - \mathbf{6} \mathbf{Q} \tag{2}$$

with 2G being the dimensionless applied pressure gradient scaled by <sup>µ</sup>*<sup>V</sup> D*<sup>2</sup> . Clearly, Q < 1 for coating flow. Consequently, −3 < G = *O*(1) < 3. The final thickness T is typically greater than half the gap for blade coating where G < 0. That the pressure gradient can be constant far upstream is easily observed in blade coating for a long flat or even angled blade. We refer the reader to the numerical studies of Iliopoulos and Scriven [17], particularly their figures for different blade angles, and Mitsoulis

and Athanasopoulos [50]. For slot coating, the pressure gradient is commonly directly related to the flow rate, which can be independent of the substrate speed [10]. In this case, G can be positive or negative, and is typically assumed to be constant in numerical studies of slot coating [15,22,25]. Finally, the value of G varies typically from −0.05 to −0.20 in blade coating [17,24,50], and from 0.05 to 0.20 in slot coating [11,18]. positive or negative, and is typically assumed to be constant in numerical studies of slot coating [15,22,25]. Finally, the value of G varies typically from −0.05 to −0.20 in blade coating [17,24,50], and from 0.05 to 0.20 in slot coating [11,18]. Figure 1 illustrates the flow configuration schematically where dimensionless notations are used. We conveniently choose the stationary wall (slot die or blade) to lie along the x axis at z = 0,

related to the flow rate, which can be independent of the substrate speed [10]. In this case, G can be

*Fluids* **2020**, *5*, x 6 of 51

flow. Consequently, 3G 1 3 *O* . The final thickness T is typically greater than half the gap for blade coating where G < 0. That the pressure gradient can be constant far upstream is easily observed in blade coating for a long flat or even angled blade. We refer the reader to the numerical

*V D* 

. Clearly, *Q* < 1 for coating

with 2*G* being the dimensionless applied pressure gradient scaled by <sup>2</sup>

Figure 1 illustrates the flow configuration schematically where dimensionless notations are used. We conveniently choose the stationary wall (slot die or blade) to lie along the x axis at z = 0, and the moving substrate at *z* = 1. The z axis lies across the channel, with the origin coinciding with the contact line (*x* = *z* = 0) between the stationary wall and the free surface z = ζ(x). The figure depicts the two possibilities for the flow far upstream that typically correspond to a thin or a thick coating film, involving an adverse (G > 0) or a favorable (G < 0) pressure gradient, respectively. In both cases, the pressure gradient is generally expected to be constant. Figure 1 also illustrates the four regions or layers that constitute the physical domain, which will be discussed in Section 2.2. and the moving substrate at *z* = 1. The z axis lies across the channel, with the origin coinciding with the contact line (*x = z* = 0) between the stationary wall and the free surface z x . The figure depicts the two possibilities for the flow far upstream that typically correspond to a thin or a thick coating film, involving an adverse G 0 or a favorable G 0 pressure gradient, respectively. In both cases, the pressure gradient is generally expected to be constant. Figure 1 also illustrates the four regions or layers that constitute the physical domain, which will be discussed in Section 2.2.

**Figure 1.** Schematic illustration of the coating flow domain. The two possible Couette–Poiseuille profiles correspond to a favorable pressure gradient <sup>1</sup> Q ,G 0 <sup>2</sup> and an adverse pressure gradient <sup>1</sup> Q ,G 0 <sup>2</sup> . The domain is subdivided in four layers: (I) the core layer, (II) the free-surface layer (x **Figure 1.** Schematic illustration of the coating flow domain. The two possible Couette–Poiseuille profiles correspond to a favorable pressure gradient Q > 1 2 , G < 0 and an adverse pressure gradient Q < 1 2 , G > 0 . The domain is subdivided in four layers: (I) the core layer, (II) the free-surface layer (x > 0), (III) the stationary wall slip layer (x < 0), and (IV) the substrate slip layer. Double arrows indicate matching directions.

#### > 0), (III) the stationary wall slip layer (x < 0), and (IV) the substrate slip layer. Double arrows indicate *2.1. Governing Equations and Boundary Conditions*

matching directions. *2.1. Governing Equations and Boundary Conditions*  The problem is formulated in terms of the stream function ψ, with u = ψ<sup>z</sup> and w = −ψ<sup>x</sup> . The stream function and velocity corresponding to the CP flow far upstream turn out to be the leading-order solution in the core layer I, and are conveniently introduced here in terms of G as

$$
\psi\_0(\mathbf{z}) = \frac{1-\mathcal{G}}{2}\mathbf{z}^2 + \frac{\mathcal{G}}{3}\mathbf{z}^3,\\
\mathbf{u}\_0(\mathbf{z}) = (1-\mathcal{G})\mathbf{z} + \mathcal{G}\mathbf{z}^2. \tag{3}
$$

 2 3 0 1G G z zz 2 3 , <sup>2</sup> <sup>0</sup> u z 1 G z Gz . (3) In the current study, the Reynolds number, Re, is assumed to be moderately large and G = O(1) at most. If Q = 1/2 (G = 0), the flow is only driven by the forward substrate translation, resulting in drag flow. For steady laminar planar flow, the non-dimensional Navier–Stokes equations take the following form:

$$\begin{aligned} \psi\_{\rm z} \psi\_{\rm xz} - \psi\_{\rm x} \psi\_{\rm zz} &= -\mathbf{p}\_{\rm x} + \frac{1}{\text{Re}} (\psi\_{\rm xxz} + \psi\_{\rm zzz}),\\ \psi\_{\rm x} \psi\_{\rm xz} - \psi\_{\rm z} \psi\_{\rm xx} &= -\mathbf{p}\_{\rm z} - \frac{1}{\text{Re}} (\psi\_{\rm xxx} + \psi\_{\rm xzz}). \end{aligned} \tag{4}$$

Except for the stress components, a subscript with respect to x or z denotes partial differentiation. For x > 0, the kinematic and dynamic boundary conditions at the free surface z = ζ(x) are

$$
\psi(\mathbf{x} > \mathbf{0}, \mathbf{z} = \zeta) = \mathbf{0}, \tag{5}
$$

$$\mathbf{p}(\mathbf{x} > \mathbf{0}, \mathbf{z} = \zeta) \zeta' - \left. \frac{1}{\text{Re}} (2\zeta'\psi\_{\text{xz}} + \psi\_{\text{xx}} - \psi\_{\text{zz}}) \right|\_{\mathbf{z} = \zeta} = -\frac{1}{\text{Re}\text{Ca}} \frac{\zeta'\zeta''}{1 + \zeta'^2} \tag{6}$$

$$\mathbf{p}(\mathbf{x} > \mathbf{0}, \mathbf{z} = \boldsymbol{\zeta}) + \frac{1}{\text{Re}} [2\psi\_{\text{xz}} + \boldsymbol{\zeta}'(\psi\_{\text{xz}} - \psi\_{\text{xx}})]\_{\mathbf{z} = \boldsymbol{\zeta}} = -\frac{1}{\text{Re}\text{Ca}} \frac{\boldsymbol{\zeta}\prime}{1 + \boldsymbol{\zeta}'^2}.\tag{7}$$

Here, a prime denotes total differentiation. In addition, the following conditions must be satisfied at the walls and far upstream:

$$
\psi\_{\mathbf{z}}(\mathbf{x}, \mathbf{z} = 1) = 1,\qquad\qquad\psi(\mathbf{x}, \mathbf{z} = 1) = 1 - \frac{\mathbf{G}}{6},\tag{8}
$$

$$
\psi\_{\mathbf{z}}(\mathbf{x} < \mathbf{0}, \mathbf{z} = \mathbf{0}) = \mathbf{0}, \tag{9}
$$

$$
\psi(\mathbf{x} < \mathbf{0}, \mathbf{z} = \mathbf{0}) = \mathbf{0}, \tag{10}
$$

$$
\psi(\mathbf{x}\rightarrow-\infty,\mathbf{z}) \sim \frac{1-\mathcal{G}}{2}\mathbf{z}^2 + \frac{\mathcal{G}}{3}\mathbf{z}^3. \tag{11}
$$

It is worth noting that in the coating literature, the pressure scale is taken as µV/D instead of the current ρ*V* 2 scale. In the former case, the pressure becomes Rep(x, z), which we shall use as such. The overall solution strategy of problem (4)–(11) along with the flow structure is summarized next.

#### *2.2. The Physical Domain and the Flow Structure*

The physical domain is shown schematically in Figure 1. As the fluid emerges out of the channel, the shear stress vanishes at the channel exit, causing *free-surface layer* II to develop along the free surface. The *core layer* I is influenced by the thin free-surface layer and the jet contraction. It is important to recall that the core layer remains predominantly of inviscid rotational character. However, this inviscid layer cannot extend to the stationary die or blade wall where significant viscous shearing occurs upstream of the stress singularity (x < 0) at the exit, forcing the formation of the *slip layer* III in the vicinity of the die or blade. Similarly, the core flow must also adjust to the shearing viscous flow in the *slip layer* IV upstream and downstream. Consequently, the predominantly inviscid core flow profile cannot satisfy adherence at the rigid walls, thus causing a slip layer to form with thickness that is not expected to be significant at high Reynolds number, but can be large otherwise. It is therefore fundamentally important to examine the flow in the slip layer and the extent of upstream influence. In fact, the thickness of the wall layers will turn out to be of O Re−<sup>1</sup> . The layer IV along the moving substrate is expected to grow rapidly with distance downstream of the channel exit.

The flow in each layer of the domain (Figure 1) is dominated by different physical mechanisms with corresponding characteristic length scales. The flow in layers II, III and IV, near the free surface, the die or blade and the moving substrate, respectively, is shear dominated, and is of the boundary-layer type. In the core layer I, both shear and elongation are in balance as a result of the predominance of the Couette–Poiseuille character of the flow upstream and downstream of the channel exit. At the channel exit, x = 0, the shear stress undergoes a step change from a non-zero value τxz(x < 0, z = 0) ≈ 1 − G at the die or blade to zero at the free surface *z* = ζ(*x*). The effect of this drop diffuses upstream inside the channel (x < 0) over a certain distance where fully developed Couette–Poiseuille flow is recovered, and downstream (x > 0), toward the moving substrate. In the limit of infinite Reynolds number, the flow retains the Couette–Poiseuille profile (3), which will turn out to be the lowest-order solution in the core layer I. In this limit, the jet height remains horizontal, with ζ(x) = 0. In this case, the jet does not contract, and mass remains conserved throughout the upstream channel and downstream jet. For finite Reynolds number, the Couette–Poiseuille profile is modified when the fluid leaves the channel in the form of the wall jet. As the fluid detaches itself from the stationary wall, the removal of the wall stress causes the free-surface layer II to develop, where the Couette–Poiseuille profile adjusts itself so as to satisfy the stress-free (zero-traction) condition at the free surface. At infinite Reynolds number, the zero-traction condition would not be imposed; there is no viscous mechanism for the stress

singularity to dissipate. In this limit, all the conditions of the problem would be satisfied since the Couette–Poiseuille profile continues undisturbed downstream from the exit. However, as the inviscid flow is governed by the Euler equations with no unique solution, the fully developed Couette–Poiseuille flow can be assumed to be everywhere the proper inviscid limit. Consequently, the core flow is, not affected by the in the free-surface layer II to lowest order; this layer, however, induces perturbations to the Couette–Poiseuille flow when higher-order terms are included, with some upstream influence. This assumption is similar to the one made for channel and tube flows with mild [51] and severe [52] constrictions, where the flow field in the core region, to leading order, satisfies the inviscid equations of motion. The reader is again referred to the book by Sobey [39] for other asymptotic analyses on interactive boundary layers.

The analysis of the flow in the core layer I is conducted separately upstream (x < 0) and downstream (x > 0) of the exit. The flow is obtained as a perturbation to the Couette–Poiseuille flow by solving Equation (4) subject to the far-upstream condition (11) and the matching condition with the free-surface layer II as the vertical double-arrow between layers I and II indicates in Figure 1. The upstream and downstream flows in layer I are matched at x = 0 as the horizontal double-arrow indicates. Boundary-layer solutions are sought in layers II, III and IV near the free surface and the walls, but not all layers admit a similarity solution. For the flow downstream of the exit, a similarity solution is possible in the free-surface layer II.

This determines completely the flow in layers I and II. The details will be given in Sections 3 and 4. We observe that no matching is required at x = 0 between the free-surface layer II and the slip layer III. Consequently, the singularity is entirely circumvented in the solution process, which constitutes a major advantage of the present approach. In Section 3, we also use a thin-film approach to analyze the flow far downstream. This is an appropriate approach since the flow is completely of the boundary-layer type at some distance downstream from the exit. In Section 5, we examine the flow near the stationary wall or slip layer III. The flow in layer III is matched onto the flow in the core layer I for x < 0 as the vertical double-arrow between layers I and III indicates. Finally, the flow in the slip layer IV near the moving substrate is examined in Section 6. The non-similarity solution is obtained separately for *x* < 0 and x > 0, and matched at x = 0.

#### **3. The Flow in the Free-Surface Layer II and the Coating Film Profile Near and Far from the Exit**

In this section, we examine the flow structure in the free-surface layer II. As the layer grows with downstream distance, it eventually invades the entire film region, at which point we adopt a thin-film approach to capture the flow far from the die or blade exit. The flows near and far downstream from the exit are matched to obtain the profile of the coating film everywhere downstream.

#### *3.1. The Flow in the Free-Surface Layer Close to the Die or Blade Exit*

In addition to the kinematic and dynamic conditions (5)–(7) at the free surface, the matching of the flow at the edge of the free-surface layer and the core layer also provides conditions needed to solve for the flow near the free surface as well as the core stream function and pressure terms. The matching process is detailed in Appendix A. For a successful application of the matching rule (A1), it is required that the stretching transformation must be expressed in the canonical form y = εη where *y* = *z* − ζ(*x*) is the near boundary transverse coordinate in the free-surface layer II. Here, ε(Re) << 1 is the small parameter in the problem or the aspect ratio as referred to by Weinstein and Ruschak [2], which will be defined precisely shortly. The core expansion in this case, must be written in terms of y, not z; otherwise (A1) will be satisfied only approximately.

To examine the free-surface layer structure, we let ε = Re−α, where α is to be determined. Anticipating that the height ζ of the free surface is of the same order of magnitude as the boundary-layer

thickness, one can write ζ(*x*) = ε*h*(*x*), and henceforth work with *h*(*x*), with the matching indicating that *h*(*x*) = O(1). The following change of coordinates is introduced, namely,

$$\mathbf{x} = \mathsf{E}, \qquad \qquad \mathbf{z} = \mathbf{y} + \mathsf{L}(\mathbf{x}) = \varepsilon(\mathfrak{n} + \mathbf{h}). \tag{12}$$

The use of ξ instead of x helps emphasize the distinction between the core region where the (x, z) variables are used, and the free-surface layer where the (ξ, η) variables are used, and makes clearer the mathematical development in each region. The aim is to find a solution of the problem in the (ξ, η) plane in the form of a boundary-layer expansion in ε. In order to match this to the core (predominantly) CP flow, it is necessary to have ψ ∼ 1−G 2 ε 2η <sup>2</sup> as <sup>η</sup> → ∞ in layer II, to lowest order in <sup>ε</sup>. Therefore, ψ must be of order ε 2 . In order to determine the value of α, the convective and viscous terms in the transformed momentum Equation (4) must balance. This is achieved upon taking the value of α = 1/3 similar to the case of a Newtonian jet [31–33] as well as a non-Newtonian jet [37]. The streamwise and transverse velocity components are now expressed in terms of the stream function as u = ψ<sup>z</sup> = <sup>1</sup> <sup>ε</sup>ψ<sup>η</sup> and w = −ψ<sup>x</sup> = −ψ<sup>ξ</sup> + h <sup>0</sup>ψη, respectively. Considering the fact that the streamwise velocity u in layer II must match the velocity in the core layer I, that is u ∼ (1 − G)εη as η → ∞, it is concluded that *u* is of order ε. From continuity, *w* is of order ε 2 . Consequently, the momentum conservation Equation (4) become

$$
\begin{split}
\psi\_{\eta}\psi\_{\xi,\eta} - \psi\_{\xi}\psi\_{\eta\eta} &= -\varepsilon^{2} \Big( \mathbf{p}\_{\xi} - \mathbf{h}' \mathbf{p}\_{\eta} \Big) \\ &+ \varepsilon^{2} \psi\_{\eta\eta\eta} + \varepsilon^{4} \Big( \psi\_{\xi,\xi,\eta} - 2\mathbf{h}' \psi\_{\xi\eta\eta} + \mathbf{h}'^{2} \psi\_{\eta\eta\eta} - \mathbf{h}'' \psi\_{\eta\eta} \Big) 
\end{split} \tag{13}
$$

$$\begin{split}-\psi\_{\eta}\psi\_{\xi,\xi} + \psi\_{\xi}\psi\_{\xi,\eta} + \mathbf{h}''\psi\_{\eta}^{2} &\quad + \mathbf{h}'\Big(\psi\_{\eta}\psi\_{\xi,\eta} - \psi\_{\xi}\psi\_{\eta\eta}\Big) \\ &= -\mathbf{p}\_{\eta} - \varepsilon^{2}\Big(\psi\_{\xi\eta\eta} - \mathbf{h}'\psi\_{\eta\eta\eta}\Big) - \varepsilon^{4}\Big(\frac{\partial}{\partial\xi} - \mathbf{h}'\frac{\partial}{\partial\eta}\Big)^{3}\psi. \end{split} \tag{14}$$

Similarly, the boundary conditions on the free surface, i.e., at η = 0, can be rewritten as

$$
\psi = 0,\tag{15}
$$

$$\mathbf{h'p} + \psi\_{\eta\eta} + \varepsilon^2 \left(\mathbf{h''}\psi\_{\eta} + \mathbf{h'^2}\psi\_{\eta\eta}\right) = -\mathbf{Ca^{-1}}\varepsilon^4 \frac{\mathbf{h'h''}}{1 + \varepsilon^2 \mathbf{h'^2}}\tag{16}$$

$$\varepsilon \mathbf{p} + \varepsilon^2 \Big( 2\psi\_{\xi\eta} - \mathbf{h}'\psi\_{\eta\eta} \Big) + \varepsilon^4 \mathbf{h}' \Big( 2\mathbf{h}'\psi\_{\eta\cdot\xi} + \mathbf{h}''\psi\_\eta - \mathbf{h}'^2\psi\_{\eta\eta} \Big) = -\mathbf{C}\mathbf{a}^{-1}\varepsilon^4 \frac{\mathbf{h}''}{1 + \varepsilon^2 \mathbf{h}'^2}.\tag{17}$$

The expansion for ψ in the free-surface layer II begins with a term in ε 2 , and precedes in powers of ε as

$$
\psi(\xi,\eta) = \varepsilon^2 \Psi\_2(\xi,\eta) + \varepsilon^3 \Psi\_3(\xi,\eta) + \cdots \ . \tag{18}
$$

Similarly, h is expanded as

h(ξ) = ε −1 ζ(ξ) = h0(ξ) + εh1(ξ) + · · · . (19)

From (18) and (19), it is concluded that the pressure p in the free-surface layer is of order ε 4 . Hence,

$$\mathbf{p}(\xi,\eta) = \varepsilon^4 \mathbf{P}\_4(\xi,\eta) + \varepsilon^5 \mathbf{P}\_5(\xi,\eta) + \cdots \,\,\,\,\,\,\,\,\tag{20}$$

This leaves the leading-order shear-stress term in Equation (16) to balance the surface tension term, indicating that surface tension is of order <sup>ε</sup> 2 Ca <sup>=</sup> <sup>1</sup> CaRe2/3 . We shall see in Section 3.2 that for a thin film at high Reynolds number, surface tension is of order <sup>1</sup> CaRe<sup>3</sup> (refer also to Weinstein and Ruschak [2]). Consequently, for coating flow, which behaves like a thin film far downstream of the exit, surface tension becomes negligible when <sup>1</sup> CaRe<sup>3</sup> <sup>&</sup>lt; 1 CaRe2/3 << <sup>1</sup> or Ca >> <sup>1</sup> Re2/3 <sup>&</sup>gt; 1 Re<sup>3</sup> . So, even for a relatively moderate Reynolds number, say Re = 64, Ca needs only be greater than 0.06 for

surface tension to be negligible. In this work, we assume a moderate surface tension with Ca remaining of order one.

Recalling that u = <sup>1</sup> <sup>ε</sup>ψ<sup>η</sup> and w = −ψ<sup>ξ</sup> + h <sup>0</sup>ψη, then the velocity components take the following for:

$$\mathfrak{u}(\xi,\mathfrak{v}) = \varepsilon \mathsf{U}\_1(\xi,\mathfrak{v}) + \varepsilon^2 \mathsf{U}\_2(\xi,\mathfrak{v}) + \cdots \mathsf{ },\tag{21}$$

$$\mathbf{w}(\xi,\mathfrak{v}) = \varepsilon^2 \mathbf{W}\_2(\xi,\mathfrak{v}) + \varepsilon^3 \mathbf{W}\_3(\xi,\mathfrak{v}) + \cdots \tag{22}$$

In this case, U<sup>1</sup> = Ψ2η, U<sup>2</sup> = Ψ3<sup>η</sup> and W<sup>2</sup> = −Ψ2<sup>ξ</sup> + h<sup>0</sup> <sup>0</sup>Ψ2η, and so on.

We next consider the leading-order solution of problem (14) and (15). Thus, to leading order in ε, Equation (13) reduce to:

$$
\Psi\_{2\eta}\Psi\_{2\xi,\eta} - \Psi\_{2\xi}\Psi\_{2\eta\eta} = \Psi\_{2\eta\eta\eta}.\tag{23}
$$

The above problem is similar to the case of symmetric free jet [31–33] with different boundary conditions. The conditions at the free surface are deduced from Equation (15) as

$$
\Psi\_2(\xi, 0) = \Psi\_{2\eta\eta}(\xi, 0) = 0. \tag{24}
$$

To complete the problem for Ψ2, another boundary condition is required, which is obtained by matching the flow at the edge of the free-surface layer with the flow in the core layer downstream of the channel exit. This is shown in Appendix A. To this order, there is no interaction between the free-surface layer II and the flow in the core layer I, this latter still retaining the Couette-Poiseuille profile. The resulting condition is conveniently written here:

$$\Psi\_2(\xi, \eta \to \infty) \sim \left(\frac{1-\mathcal{G}}{2}\right)\eta^2. \tag{25}$$

We observe that the growth of the free-surface layer II can be estimated by considering the balance of the viscous and inertial terms in Equation (23) and using Equation (25), thus indicating that η ∝ ξ Ψ<sup>2</sup> ∼ 2 1−G ξ η2 . This leads to y ∝ 2 1−G 1/3 x Re 1/3 or

$$
\delta \propto \left(\frac{2}{1-\mathcal{G}} \frac{\chi}{\text{Re}}\right)^{1/3},\tag{26}
$$

which suggests that the thickness of the free-surface layer for a coating film decreases like Re−1/3 with increasing Reynolds number. However, this also suggests that, for the same substrate speed (same Reynolds number), the free-surface layer is expected to be thicker for slot coating (G > 0) compared to blade coating (G < 0). As an interesting consequence, the surface layer grows faster for the thinner film Q < 1 2 in slot coating than for the thicker film Q > 1 2 in blade coating. This fundamental observation has an important practical significance in modelling coating film flow, which will be discussed in Section 3.2 when the near exit flow is matched to the thin-film flow. Finally, as we shall see, the boundary-layer growth is different along the stationary and moving walls.

We now turn to the solution of problem (23)–(25), which admits a similarity profile: Ψ2(ξ, η) = ξ 2/3 f2(θ), where θ = ηξ−1/3 is the similarity variable. The problem for f2(θ) is given by

$$\text{'}\text{'}\text{'}\text{''}\_2 + \text{2f}\_2\text{'}\text{'}-\text{'}\_2^2 = 0,\tag{27}$$

$$\mathbf{f\_2(0)} = \mathbf{f\_2''(0)} = \mathbf{0},\tag{28}$$

$$\text{tf}\_2(\theta \to \infty) \sim \left(\frac{1-\mathcal{G}}{2}\right)(\theta + \mathbf{c}\_1)^2,\tag{29}$$

where c<sup>1</sup> = (1 − G) −1/3 d<sup>1</sup> is a constant function of G to be determined numerically. Problem (27)–(29) is solved as an initial-value problem using a fourth-order Runge–Kutta scheme (IMSL-DIVERK), coupled with a shooting technique. Equation (27) is integrated subject to conditions (28) and a guessed

value f 0 2 (0) for the slope at the origin. The slope is adjusted until reasonable matching is achieved with the asymptotic form (29) at large θ, or, more precisely, until f 00 2 ∼ 1 − G is reached. The value of c<sup>1</sup> is then determined upon matching the numerical solution and its asymptotic form. Matching the numerical solution with its asymptotic behavior, we find d<sup>1</sup> = 0.892. Moreover, the slope f 0 2 (0) = (1 − G) 2/3 e1 is related to the velocity at the film surface where we numerically find e<sup>1</sup> = 1.611. We observe that the numerical solution was carried out by scaling out the factor 1–G from the problem. The variables were rescaled by letting f2(θ) = (1 − G) 1/3 g2 (t) and t = (1 − G) 1/3 θ. In this case, problem (27)–(29) reduces to: 3g<sup>000</sup> 2 + 2g<sup>2</sup> g 00 2 − g 02 <sup>2</sup> = 0, g<sup>2</sup> (0) = g 00 2 (0) = 0, g<sup>2</sup> (t → ∞) ∼ (t+d<sup>1</sup> ) 2 2 . Pursuing the solution and the matching process to next order, the governing equation for Ψ3(ξ, η) becomes

$$\mathbf{V}\_{2\eta}\mathbf{V}\_{3\xi\eta} + \mathbf{V}\_{3\eta}\mathbf{V}\_{2\xi\eta} - \mathbf{V}\_{2\xi}\mathbf{V}\_{3\eta\eta} - \mathbf{V}\_{3\xi}\mathbf{V}\_{2\eta\eta} = \mathbf{V}\_{3\eta\eta\eta} \tag{30}$$

subject to the boundary conditions

$$
\Psi\_{\mathcal{B}}(\xi, 0) = \Psi\_{\mathcal{B}^{\eta \eta \eta}}(\xi, 0) = 0,\tag{31}
$$

$$\Psi\_{\mathcal{B}}(\xi,\mathfrak{u}) \sim \frac{G}{3}\mathfrak{n}^3 \text{ as } \mathfrak{n} \to \infty. \tag{32}$$

A similarity solution is also possible, namely Ψ3(ξ, η) = ξ f3(θ), resulting in the following problem for f3(θ):

$$\mathfrak{H}\_{3}^{\prime\prime\prime} + \mathfrak{H}\_{2}\mathfrak{f}\_{3}^{\prime\prime} - \mathfrak{H}\_{2}^{\prime}\mathfrak{f}\_{3}^{\prime} + \mathfrak{H}\_{2}^{\prime\prime}\mathfrak{f}\_{3} = 0,\tag{33}$$

$$\mathbf{f\_3(0)} = \mathbf{f\_3''(0)} = \mathbf{0},\tag{34}$$

$$\text{rf}\_3(\theta \to \infty) \sim \frac{\text{G}}{3} \Big[ \left( \theta + \mathbf{c}\_1 \right)^3 - \frac{6}{1 - \mathbf{G}} \Big] + \mathbf{c}\_2(\theta + \mathbf{c}\_1). \tag{35}$$

here c<sup>2</sup> = G(1 − G) −2/3 d2, where the constant d<sup>2</sup> = 1.654 is determined numerically. Here again, the numerical solution is obtained for the rescaled equation and boundary conditions:

$$2\mathbf{g\_3'''} + 2\mathbf{g\_2g\_3''} - 3\mathbf{g\_2g\_3'} + 3\mathbf{g\_2''}\mathbf{g\_3} = \mathbf{0},\tag{36}$$

$$\mathbf{g\_3(0) = g\_3""(0) = 0, \ g\_3(t \to \infty) \sim \frac{1}{3} \left[ (t + \mathbf{d\_1})^3 - 6 \right] + \mathbf{d\_2(t + \mathbf{d\_1})}.\tag{37}$$

We find it is helpful to summarize in Table 1 the constants and numerical values that are used in the numerical results.

**Table 1.** Constants arising in the solution of f2(θ) and f3(θ).


We recall that G is directly related to the flow rate (or coating thickness) through G = 3(1 − 2Q). The height of the free surface is then obtained to the current order by substituting (A11) and (A18) into Equation (19), yielding

$$\zeta(\mathbf{x}) = \varepsilon \mathbf{c}\_1 \mathbf{x}^{1/3} + \varepsilon^2 \frac{\mathbf{c}\_2}{1 - \mathbf{G}} \mathbf{x}^{2/3} = \mathbf{c}\_1 \left(\frac{\mathbf{x}}{\text{Re}}\right)^{1/3} + \frac{\mathbf{c}\_2}{1 - \mathbf{G}} \left(\frac{\mathbf{x}}{\text{Re}}\right)^{2/3}. \tag{38}$$

This expression reveals clearly the intricate interplay between inertia and flow rate (coating thickness). In particular, the free-surface height ζ(x) or the thickness χ(x) ≡ 1 − ζ(x) can exhibit an extremum depending on the value and signs of the constants involved. An important observation to make here is the absence of surface-tension effect on the film profile given by Equation (38). This is not surprising since Ca = *O*(1), which is relatively large. In this range of capillary number, the film thickness remains essentially independent of surface tension as this has been clearly demonstrated experimentally. The reader is referred, among others, to Lee et al. [6] and Becker and Wang [7] where film thickness measurements were reported for slot coating. In their numerical simulation.

Saito and Scriven [15] examined the influence of flow parameters, such as the flow rate, the capillary number and the Reynolds number, on the shape of the meniscus in slot coating. They observed that the effect of inertia was the least evident to interpret given the non-monotonic response in the meniscus profiles as they varied the Reynolds number. A similar non-monotonic response was also reported more recently by Carvalho and Kheshgi [9] who plotted the numerically found film profiles for different capillary numbers for constant Re/Ca ratio. Incidentally, this ratio is the Ohnesorge number Oh−<sup>2</sup> <sup>≡</sup> Re Ca [53], and not the "Property" number as sometime referred to in the coating literature. The non-monotonic response becomes particularly clear when the separation angle and the radius of curvature near the contact line are plotted against Re, displaying a maximum and a minimum, respectively, which is illustrated in Saito and Scriven [15].

At the free surface (*z* = ζ), the velocity becomes

$$\mathbf{u}(\mathbf{x}, \mathbf{z} = \boldsymbol{\zeta}) = \varepsilon \mathbf{x}^{1/3} \mathbf{f}\_2'(0) + \varepsilon^2 \mathbf{x}^{2/3} \mathbf{f}\_3'(0) = \left(\frac{\mathbf{x}}{\mathrm{Re}}\right)^{1/3} \mathbf{f}\_2'(0) + \left(\frac{\mathbf{x}}{\mathrm{Re}}\right)^{2/3} \mathbf{f}\_3'(0). \tag{39}$$

We note that f0 3 (0) = G(1 − G) −2/3 e<sup>2</sup> where e<sup>2</sup> = 2.85. Both the initial slope f<sup>0</sup> 3 (0) and c<sup>2</sup> increase with positive G values, reflecting a higher order strengthening effect of the adverse pressure gradient on the flow near the free surface. In contrast, the trend is reversed for negative G, pointing to a higher order weakening effect of the favorable pressure gradient on the flow close to the free surface.

It is interesting to observe that both the free-surface height (or film thickness) and the velocity depend on <sup>x</sup> Re rather than x. This is also the case for the pressure along the free surface of the meniscus, which is determined by inserting Equations (18)–(20) into condition (17). To the current order, the pressure may be written as

$$\mathrm{Re}^2 \mathbf{p}(\xi, 0) = -\frac{2}{3} \left[ \left( \frac{\mathbf{x}}{\mathrm{Re}} \right)^{-2/3} \mathbf{f}\_2'(0) + 2 \left( \frac{\mathbf{x}}{\mathrm{Re}} \right)^{-1/3} \mathbf{f}\_3'(0) \right] + \frac{2}{9} \mathrm{We}^{-1} \left( \mathrm{c}\_l \left( \frac{\mathbf{x}}{\mathrm{Re}} \right)^{-5/3} + \frac{\mathbf{c}\_2}{1 - \mathbf{G}} \left( \frac{\mathbf{x}}{\mathrm{Re}} \right)^{-4/3} \right). \tag{40}$$

Here We = Ca/Re is the Weber number. Thus, near the exit, the pressure behaves like x −2/3 in the absence of surface tension, which is different from the x <sup>−</sup>1/2 behavior reported by Aidun [14]. Moreover, the strength of the pressure singularity depends intricately on inertia, surface tension and the flow rate (or upstream pressure gradient). In the absence of surface tension, the singularity weakens with inertia and the flow rate for both slot and blade coating. As mentioned earlier. Chang et al. [8] observed that viscous and surface-tension effects become negligible for Re > 20. One expects the pressure to be vanishingly small along the free surface. In this case, Equation (40) gives an estimate of pressure on the order of Re−4/3 .

Although surface tension can play a significant role, not only by becoming dominant near the exit, but also by altering the force behavior as we shall see when we discuss the normal stress distribution shortly. For now, we consider the effects of inertia and flow rate on the film height, the velocity and the pressure along the film surface.

Figure 2 depicts the dependence of the film height (Figure 2a), the velocity (Figure 2b) and the pressure (Figure 2c) along the free surface, on the flow rate (or pressure gradient) and inertia. It is evident from the figure that as the dimensionless flow rate Q increases, the film height in Figure 2a drops as expected, signalling a thicker film and a longer relaxation distance. This behavior agrees well with existing numerical predictions as confirmed by comparing with Saito and Scriven [15] and Carvalho and Kheshgi [9]. Very close to the exit, the film height (Figure 2a), the velocity (Figure 2b) as well as the pressure (Figure 2c) increase sharply and monotonically with distance, for any flow rate, indicating a strong contraction. This behavior is the result of the dominance of the x 1/3 term in Equations (38) and (39) near the origin. Although the height and the streamwise velocity components

are continuous at the exit (as they remain equal to zero immediately before x = 0), the flow behavior is intrinsically singular through the surface slope, transverse velocity and stress components, and is well illustrated in the pressure behavior near x = 0 (Figure 2c). We observe from Figure 2a that the separation angle is consistently 90 degrees for the order of capillary and Reynolds numbers considered here. The separation angle is the angle between the normal to the stationary slot or blade wall and the normal to the surface at the contact line (x = z = 0). Further downstream, only the film at low flow rate, typically depicting the behavior of a film in slot coating, continues to grow, resulting in a smaller coating thickness. While the film height grows monotonically with distance for a drag film (Q = 0.5), the free-surface height and velocity for blade coating (Q > 0.5) exhibit a maximum before decaying, indicating a local expansion. However, the presence of the maximum should be interpreted with some caution. The non-monotonic behavior suggests a relatively strong influence of the higher-order x 2/3 term, which should not dominate the leading-order x 1/3 term if the expansions (38) and (39) are to remain uniformly valid. However, this does not seem to be the case for the range of flow rates reported in Figure 2. The expansion displayed for the thicker film in blade coating (Q > 0.5) appears to be physically real as can be established from the numerical values of the various terms. However, the validity of (38) and (39) also depends on the thickness of the free-surface layer II (boundary layer along the free surface), which is assumed to be small compared to the free-surface height. The curves in Figure 2a are based on the near-exit solution (38). The error is O x Re 2/3 , except for the Q = 0.5 (G = 0), which is O x Re 1/3 since c2(G = 0) = 0 as per Table 1. Despite the relatively large error, the case Q = 0.5 is included to illustrate the trend as Q increases. The more accurate prediction is given in Section 3.2 when the far field is examined. The pressure profiles in Figure 2c are essentially monotonic; the pressure increases and decays asymptotically to zero. A slight overshoot ahead of the decay is observed for the thicker coating films. Despite the clear dependence of the elongation rate, the second term on the left of (17) that yields the wild non-monotonic dynamics in Figure 2b, the monotonicity of the pressure is mainly caused by surface tension. Although this behavior is reported for We = 1, the same trend is predicted essentially for any moderate level of surface tension. This behavior sharply contrasts that of the normal stress, which depends strongly on surface tension as we shall now see. *Fluids* **2020**, *5*, x 15 of 51 various terms. However, the validity of (38) and (39) also depends on the thickness of the free-surface layer II (boundary layer along the free surface), which is assumed to be small compared to the freesurface height. The curves in Figure 2a are based on the near-exit solution (38). The error is 2/3 <sup>x</sup> <sup>O</sup> Re , except for the Q = 0.5 (G = 0), which is 1/3 <sup>x</sup> <sup>O</sup> Re since cG0 0 <sup>2</sup> as per Table 1. Despite the relatively large error, the case Q = 0.5 is included to illustrate the trend as Q increases. The more accurate prediction is given in Section 3.2 when the far field is examined. The pressure profiles in Figure 2c are essentially monotonic; the pressure increases and decays asymptotically to zero. A slight overshoot ahead of the decay is observed for the thicker coating films. Despite the clear dependence of the elongation rate, the second term on the left of (17) that yields the wild non-monotonic dynamics in Figure 2b, the monotonicity of the pressure is mainly caused by surface tension. Although this behavior is reported for We = 1, the same trend is predicted essentially for any moderate level of surface tension. This behavior sharply contrasts that of the normal stress, which depends strongly on surface tension as we shall now see.

**Figure 2.** *Cont.*

*Fluids* **2020**, *5*, x 16 of 51

**Figure 2.** Meniscus profiles (**a**), corresponding streamwise velocity (**b**) and pressure (**c**) along the free surface for various coating thickness or flow rates for We = 1. The drag-out case (Q = 0.5) is also included for reference. **Figure 2.** Meniscus profiles (**a**), corresponding streamwise velocity (**b**) and pressure (**c**) along the free surface for various coating thickness or flow rates for We = 1. The drag-out case (Q = 0.5) is also included for reference.

The influence of flow rate and surface tension on the (negative) primary normal stress difference tt nn (Hoop stress) along the free surface is displayed in Figure 3 for Re = 10. Both the cases of slot coating (Q < 0.5) and blade coating (Q > 0.5) are illustrated in the presence (Figure 3a) and the absence (Figure 3b) of surface tension. Intricate non-monotonic behavior can be expected for blade coating, resulting from the dynamics already observed in Figure 2a,b, thus allowing the possibility of extrema, and the waviness displayed in the free surface and velocity profiles. Indeed, Figure 3a shows that the normal stress exhibits a relatively strong localized maximum followed by an asymptotic decay towards zero further downstream. This behavior is reminiscent of that reported by Bajaj et al. [19] for viscoelastic creeping flow. They examined numerically the influence of The influence of flow rate and surface tension on the (negative) primary normal stress difference τtt − τnn (Hoop stress) along the free surface is displayed in Figure 3 for Re = 10. Both the cases of slot coating (Q < 0.5) and blade coating (Q > 0.5) are illustrated in the presence (Figure 3a) and the absence (Figure 3b) of surface tension. Intricate non-monotonic behavior can be expected for blade coating, resulting from the dynamics already observed in Figure 2a,b, thus allowing the possibility of extrema, and the waviness displayed in the free surface and velocity profiles. Indeed, Figure 3a shows that the normal stress exhibits a relatively strong localized maximum followed by an asymptotic decay towards zero further downstream. This behavior is reminiscent of that reported by Bajaj et al. [19] for viscoelastic creeping flow. They examined numerically the influence of Weissenberg number and viscosity ratio on the normal stress along the free surface for dilute polymer solutions. Their figures

also show a local maximum for τtt − τnn close to the die exit followed by a decay similar to Figure 3a. A major difference, however, is worth noting here: the plots of Bajaj et al. [19] do not clearly display the singularity at the origin. They do assess, on the other hand, the existence and strength of the singularity by examining the behavior of the rate-of-strain components. As to the crucial role of surface tension, especially near the exit, it is demonstrated in Figure 3b, where the Hoop stress is plotted in the absence of surface tension. Comparison between Figure 3a,b clearly indicates that the dynamics exhibited in the stress near the origin in Figure 3a are caused by surface tension. Figure 3b shows that the Hoop stress increases essentially monotonically, exhibiting a very weak overshoot for the thicker coating film before decaying to zero. Weissenberg number and viscosity ratio on the normal stress along the free surface for dilute polymer solutions. Their figures also show a local maximum for tt nn close to the die exit followed by a decay similar to Figure 3a. A major difference, however, is worth noting here: the plots of Bajaj et al. [19] do not clearly display the singularity at the origin. They do assess, on the other hand, the existence and strength of the singularity by examining the behavior of the rate-of-strain components. As to the crucial role of surface tension, especially near the exit, it is demonstrated in Figure 3b, where the Hoop stress is plotted in the absence of surface tension. Comparison between Figure 3a,b clearly indicates that the dynamics exhibited in the stress near the origin in Figure 3a are caused by surface tension. Figure 3b shows that the Hoop stress increases essentially monotonically, exhibiting a very weak overshoot for the thicker coating film before decaying to zero.

**Figure 3.** Influence of flow rate on the normal-stress difference along the film surface as function of distance at Re = 10, and (**a**) We = 5 and (**b**) in the absence of surface tension.

Figure 4 displays the variation with the inclination angle for Re =10 and We = 0.1 of the free-surface curvature (Figure 4a) and pressure along the free surface (Figure 4b). The abscissa is the angle of inclination of the normal to the free surface from the z direction, as depicted in the inset of Figure 4b. We see that the curvature rises with increasing flow rate. This is expected as the meniscus in slot coating contracts more than in blade coating. The results in Figure 4 are overall qualitatively similar to the numerical results of Saito and Scriven [15] and Lee et al. [21], when compared to their figures for the slot coating flow. The agreement is particularly obvious when comparison is made with the high Ca curves of Saito and Scriven [15]. Similar to their curve for Ca = 2, we also find that the surface curvature exhibits a change in concavity at low inclination angle for any flow rate. At low capillary number, Saito and Scriven [15] as well as Lee et al. [21] predict a linear growth, which remains monotonic for very low capillary number. Saito and Scriven observed that a large portion of the curvature is constant when the capillary number is very small; but the curvature varies more and more rapidly, displaying a maximum, as Ca increases, a behavior that is also reflected in the rapid variation and the maxima we see in Figure 4a, especially for blade coating (Q > 0.4). Simultaneously, the pressure in Figure 4b displays a similar linear growth for small inclination angle as in the earlier numerical studies, but tends to increase sharply for larger angles as the singularity is approached near φ = 90 degrees. This behavior is also captured by Lee et al. [21].

The similarity between our high-Re results and the high-Ca results at Re = 0 of Saito and Scriven [15] highlights the crucial role of surface tension. In this regard, it is helpful to examine the current predictions for the curvature relative to a flow with higher surface tension. For Ca < 1, the shape of the meniscus near the separation line has often been approximated as a static meniscus, which was deduced by Ruschack (1976) using a quasi-static approach where the shape of the meniscus is the arc of a circle of constant curvature as predicted by Landau and Levich (1942). The theory imposes an upper bound on the meniscus curvature as it cannot exceed <sup>2</sup> 1−Q (when using our dimensionless notations). Based on the results in Figure 4a, this criterion appears to be plausible at best for the lowest flow rate considered here (Q = 0.4). This observation corroborates well that of Saito and Scriven [15] who found that the quasi-static assumption is valid only for large surface tension, which is not the case here. The convective (dynamic) effects are simply too dominant for Ruschak's approximate to approximation hold.

#### *3.2. The Coating Profile near and Far from the Exit*

In order to capture the coating profile at any location downstream of the exit, we exploit the flow structure as the free-surface layer II grows and invades the entire film region, at a critical location x = xc. Downstream of this location, the flow becomes of the boundary-layer or thin-film type. We exploit this simplification and formulate the flow using a thin-film approach, and match it at x = x<sup>c</sup> with the flow obtained earlier near the exit. Far from the exit for x > xc, the conservation equations, adherence and no-penetration conditions at the moving substrate z = 1, and the kinematic and dynamic conditions at the free surface z = ζ(x) reduce to

$$\mathbf{u}\_{\mathbf{x}} + \mathbf{w}\_{\mathbf{z}} = \mathbf{0},\\\mathrm{Re}(\mathbf{u}\mathbf{u}\_{\mathbf{x}} + \mathbf{w}\mathbf{u}\_{\mathbf{z}}) = -\mathrm{Re}\mathbf{p}\_{\mathbf{x}} + \mathbf{u}\_{\mathbf{z}\mathbf{z}\prime}\mathbf{p}\_{\mathbf{z}} = \mathbf{0},\tag{41}$$

$$\mathbf{u}(\mathbf{x}, \mathbf{z} = 1) = 1, \; \mathbf{w}(\mathbf{x}, \mathbf{z} = 1) = 0,\tag{42}$$

$$\mathbf{w}(\mathbf{x}, \mathsf{L}) = \mathbf{u}(\mathbf{x}, \mathsf{L}) \mathsf{L}'(\mathbf{x}), \tag{43}$$

$$\mathfrak{u}\_{\mathbf{z}}(\mathbf{x},\zeta) = \mathfrak{0} \ \mathfrak{p}(\mathbf{x},\zeta) = \frac{\zeta'(\mathbf{x})}{\mathrm{Re}\mathbf{C}\mathbf{a}}.\tag{44}$$

*Fluids* **2020**, *5*, 180

We proceed by using a Karman–Pohlhausen approach, and integrate the momentum equation across the film, using condition (43) and noting that p<sup>x</sup> (x) = ζ 00 (x) ReCa , to obtain

*Fluids* **2020**, *5*, x 19 of 51

$$\frac{d}{d\mathbf{x}}\int\_{\zeta}^{1} \mathbf{u}^{2} d\mathbf{z} = -\frac{\zeta''(1-\zeta)}{\mathrm{Re}\mathbf{C}\mathbf{a}} + \mathrm{Re}^{-1}\mathbf{u}\_{\mathbf{z}}(\mathbf{x}, \mathbf{z} = 1). \tag{45}$$

**Figure 4.** Variation of the meniscus curvature (**a**) and the free-surface pressure (**b**) with the freesurface inclination angle for various flow rates for slot coating Q 0.5 and blade coating **Figure 4.** Variation of the meniscus curvature (**a**) and the free-surface pressure (**b**) with the free-surface inclination angle for various flow rates for slot coating (Q ≤ 0.5) and blade coating (Q ≥ 0.5), for Re = 10 and We = 0.1.

In order to capture the coating profile at any location downstream of the exit, we exploit the flow structure as the free-surface layer II grows and invades the entire film region, at a critical location <sup>c</sup> x x . Downstream of this location, the flow becomes of the boundary-layer or thin-film type. We exploit this simplification and formulate the flow using a thin-film approach, and match it at <sup>c</sup> x x

Q 0.5 , for Re = 10 and We 0.1 .

*3.2. The Coating Profile Near and Far from the Exit* 

*Fluids* **2020**, *5*, 180

Different levels of accuracy have been adopted in the literature for the velocity profile across the film in the presence of inertia, ranging from the simplest parabolic, cubic and quartic profiles [2] to spectral expansions [54,55]. Another alternative would be the use of the asymptotic approach of Higgins and Scriven [40] for visco-capillary flow, which is based on a small departure from the plug flow that prevails far downstream. For our purpose here, we adopt the parabolic profile, which is the simplest and most commonly used, to obtain the following velocity distribution in terms of the thickness χ(x) ≡ 1 − ζ(x) as

$$\mathbf{u}(\mathbf{x},\mathbf{z}) = -\frac{3}{2} \frac{\chi\_{\rm os} - \chi(\mathbf{x})}{\chi^3(\mathbf{x})} \left[ (1 - \mathbf{z})^2 - 2\chi(\mathbf{x})(1 - \mathbf{z}) \right] + 1,\tag{46}$$

where χ<sup>∞</sup> = χ(x → ∞) = Q is the coating film thickness far downstream. The velocity along the free surface is then u(x, z = ζ) = <sup>1</sup> 2 3 χ<sup>∞</sup> χ − 1 . Substituting (45) into Equation (46) leads to the following equation for the thickness of the coating meniscus:

$$\frac{5}{\text{Re}\,\text{Ca}}\chi^3\frac{\text{d}^2\chi}{\text{d}\chi^2} - \left(\chi^2 - 6\chi\_\infty^2\right)\frac{\text{d}\chi}{\text{d}\chi} + \frac{15}{\text{Re}}(\chi - \chi\_\infty) = 0. \tag{47}$$

This equation requires two boundary conditions, which can be imposed as initial conditions at the location x = xc, where the near-exit profile given by (38) matches the far-exit profile dictated by (47). For this, we choose to match the surface height and its slope. This leaves x<sup>c</sup> as unknown, which we determine by matching the concavity. In the absence of surface tension, only the height and slope need to be matched. In this case, (47) admits an analytical solution:

$$\chi^2 - \chi\_{\rm c}^2 + 2\chi\_{\rm co}(\chi - \chi\_{\rm c}) - 10\chi\_{\rm co}^2 \ln\left(\frac{\chi - \chi\_{\rm co}}{\chi\_{\rm c} - \chi\_{\rm co}}\right) = \frac{\mathfrak{A}0}{\mathrm{Re}}(\chi - \chi\_{\rm c}),\tag{48}$$

where χ<sup>c</sup> ≡ χ(x = xc) is the thickness at xc. In this case, upon matching the thickness and the slope, x<sup>c</sup> and χ<sup>c</sup> are determined by solving the following two equations:

$$\chi\_{\rm c} = 1 - \mathbf{c}\_{\rm l} \left(\frac{\chi\_{\rm c}}{\rm Re}\right)^{1/3} - \frac{\mathbf{c}\_{\rm 2}}{(1 - \mathbf{G})} \left(\frac{\chi\_{\rm c}}{\rm Re}\right)^{2/3} \tag{49}$$

$$45\frac{\chi\_{\rm c} - \chi\_{\rm co}}{\chi\_{\rm c}^2 - 6\chi\_{\rm co}^2} + c\_1 \left(\frac{\chi\_{\rm c}}{\rm Re}\right)^{-2/3} + 2\frac{c\_2}{(1 - G)} \left(\frac{\chi\_{\rm c}}{\rm Re}\right)^{-1/3} = 0. \tag{50}$$

It is not difficult to deduce that the matching location turns out to be simply proportional to Re, whereas the corresponding thickness is independent of Re. Both quantities are of course functions of the flow rate (minimum coating thickness).

Figure 5 depicts the composite thickness (Figure 5a) and the free-surface height (Figure 5b) profiles with downstream position for a typical minimum blade coating thickness range (0.5 ≤ Q = χ<sup>∞</sup> ≤ 0.9) in the absence of surface tension. Figure 5a depicts the typical monotonic decrease in thickness with distance. The profiles appear to saturate as Q approaches one. The matching location (see Table 2) x<sup>c</sup> for each profile coincides with the intersection of the free-surface layer with z = 1, which suggests that matching occurs further downstream for a thicker film, which is reflected in Figure 5b where the free-surface (boundary) layer is plotted along with the surface height. Again, the matching location coincides with the free-surface layer II reaching the moving substrate (z = 1). We observe that G drops from 0 to −2.4 when Q increases from 0.5 to 0.9, indicating that a drop in pressure in the blade region yields a thicker coating film and thinner free-surface layer. The thinning of the free-surface layer with (negatively) increasing pressure gradient is expected as the favorable pressure difference effectively contributes additional flow inertia in blade coating, as illustrated by Iliopoulos and Scriven [17].


*Fluids* **2020**, *5*, x 22 of 51

**Table 2.** Influence of flow rate on matching location and corresponding film thickness between the flow near the exit and the flow far from the exit.

**Figure 5.** Composite thickness profiles for near and far distance downstream of the exit (**a**), and corresponding free-surface heights close to the exit (**b**) for various minimum coating thickness Q 0.5 . Also shown in (**b**) are the corresponding free-surface layer heights in lighter plots. **Figure 5.** Composite thickness profiles for near and far distance downstream of the exit (**a**), and corresponding free-surface heights close to the exit (**b**) for various minimum coating thickness (Q ≥ 0.5). Also shown in (**b**) are the corresponding free-surface layer heights in lighter plots.

 xz c <sup>2</sup>

x x ,z 1 3

is generally deduced from

examining the shear stress at the moving substrate. The shear stress along the wall far downstream

x

x

. (51)

*Fluids* **2020**, *5*, 180

An estimate of the error resulting from the use of the parabolic profile can be obtained by examining the shear stress at the moving substrate. The shear stress along the wall far downstream is generally deduced from

$$\pi\_{\mathbf{x}\mathbf{z}}(\mathbf{x}\geq\mathbf{x}\_{\mathbf{c}},\mathbf{z}=1) = \mathfrak{Z}\frac{\chi(\mathbf{x})-\chi\_{\mathbf{c}\mathbf{s}}}{\chi^{2}(\mathbf{x})}.\tag{51}$$

We consider the case of drag flow, in the absence of pressure gradient. In this case, τxz(x, z = 1) = <sup>3</sup> 2 2χ(x)−1 χ2(x) . Matching the shear stress yields the following equation for the film thickness at x = xc: 2χ 2 <sup>c</sup> − 6χ<sup>c</sup> + 3 = 0, leading to the only admissible root: χ<sup>c</sup> = 3− √ 3 <sup>2</sup> = 0.634, which is close to 0.769 (from Table 2) but not exactly the same as a result of the parabolic approximation (46), suggesting an error of 20%. Consequently, the thin-film Equation (47) is based on a crude approximation of the parabolic velocity profile (see discussion by Weinstein and Ruschak [2]). Higher-order polynomials or spectral representation may be used for a more accurate description [54,55]. Our aim here, however, is to demonstrate how the current asymptotic theory, which is valid upstream of the exit and downstream close the exit, can be matched at some location with the thin-film flow to provide a formulation to predict the flow anywhere. A more thorough approach involves adopting the asymptotic flow as initial condition for a computational (finite-element) implementation, thus avoiding having to deal numerically with the singularity at the exit.

#### **4. The Flow in the Core Layer I**

The flow in the core layer (region I in Figure 1) remains predominantly of inviscid rotational character. For this reason, the flow retains predominantly the CP profile since there is little viscous mechanism for it to develop as it approaches the exit and moves beyond the exit. The core layer is central in the current formulation as it is matched at the edge of each of the three boundary layers: the free-surface layer (region II), the stationary wall layer near a slot die or blade (region III), and the moving-wall layer near the substrate or web (region IV). Since the treatment is similar to earlier studies on the symmetric jet, only a summary of the formulation is outlined here. In fact, the present core formulation reduces to that of Tillett [31] by setting G = −1. The formulation of Khayat [37] is also recovered in the limit of a Newtonian fluid. Since the core flow does not satisfy adherence, the core formulations for the current wall jet and the symmetric jet are not affected by the nature of the boundary at z = 1, irrespective of whether it is a solid line (for the wall jet) or a symmetry line (for the symmetric jet).

The core layer is divided into two different regions: upstream (x < 0) and downstream (x > 0) of the exit at x = 0. The solution of problems (4) is sought in each region separately. As in the case of the symmetric free jet [31,37], we determine ψ and p by considering a correction to the base Couette-Poiseuille flow resulting from its interaction with the free-surface layer (II). In this case, the dominant corrections for both the stream function and the pressure turn out to be of O Re−<sup>1</sup> so that

$$
\psi(\mathbf{x}, \mathbf{z}) = \psi\_0(\mathbf{x}, \mathbf{z}) + \mathrm{Re}^{-1}\overline{\psi}(\mathbf{x}, \mathbf{z}), \\
\mathbf{p}(\mathbf{x}, \mathbf{z}) = \mathrm{Re}^{-1}\overline{\mathbf{p}}(\mathbf{x}, \mathbf{z}).\tag{52}
$$

Here, we recall that ψ<sup>0</sup> = 1−G 2 z <sup>2</sup> + <sup>G</sup> 3 z 3 is the Couette-Poiseuille stream function given in (1). In this case, ψ(x, z) represents the deviation from the base flow due to its interaction between the core and free-surface layer. The character of the base flow is similar to the laminar flow of a free jet [31] and channel or tube flow with constriction [51,52,56] at high Reynolds number. In such cases, as well, the fully developed (Poiseuille) profile is the flow to leading order. When examining the flow with severe constriction, Smith [56] expanded (1) as in (52) above. Although Smith's leading order terms in (1) still satisfy the inviscid equations of motion, they do not exactly correspond to fully developed flow as in (52). Smith adopted the free-streamline theory, and obtained the proper inviscid limiting form of the Navier–Stokes equations [56].

We observe that, given the ellipticity of the governing equations in the streamwise direction, the departure from the Couette–Poiseuille profile extends upstream to include the region x < 0, causing the upstream influence. Substituting (52) into (4) leads to

$$
\mathbf{u}\_0 \overline{\boldsymbol{\psi}}\_{\mathbf{x}\mathbf{z}} - \mathbf{u}\_{0\mathbf{z}} \overline{\boldsymbol{\psi}}\_{\mathbf{x}} = -\overline{\mathbf{p}}\_{\mathbf{x}} + 2\mathbf{G} \text{ } \mathbf{u}\_0 \overline{\boldsymbol{\psi}}\_{\mathbf{x}\mathbf{x}} = \overline{\mathbf{p}}\_{\mathbf{z}}.\tag{53}
$$

Upon eliminating the pressure and noting that w = −ψ<sup>x</sup> , one deduces the following problem for w(x, z) over the ranges −∞ < x < ∞ and 0 ≤ z ≤ 1:

$$
\overline{\mathbf{w}}\_{\text{xx}} + \overline{\mathbf{w}}\_{\text{zz}} - \frac{\mathbf{u}\_0''}{\mathbf{u}\_0} \overline{\mathbf{w}} = \mathbf{0},
\tag{54}
$$

$$
\overline{\mathbf{w}}(\mathbf{x}\rightarrow -\infty, \mathbf{z}) \sim \mathbf{0},
\tag{55}
$$

$$
\overline{\mathbf{w}}(\mathbf{x}<0,\mathbf{z}=0) = \overline{\mathbf{w}}(-\infty < \mathbf{x} < \infty, \mathbf{z}=1) = \mathbf{0},\tag{56}
$$

$$\sqrt{\mathbf{w}}(\mathbf{x} > \mathbf{0}, \mathbf{z} \to \mathbf{0}) \sim \frac{2\mathbf{G}}{1 - \mathbf{G}}.\tag{57}$$

In addition, w must remain bounded as |x| → ∞. Condition (57) is deduced from (A22), which is obtained from the matching between the core and the free-surface layers. The solution of problem (54)–(57) will now be examined separately for x < 0 and x > 0.

For x < 0, it is not difficult to see that problem (54)–(57) admits a general solution of the following form [31,37]:

$$\overline{\mathbf{w}}(\mathbf{x}<0,\mathbf{z}) = -\sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{e}^{\beta\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}(\mathbf{z}).\tag{58}$$

Here, Vn(z) and β<sup>n</sup> are orthogonal shape functions and corresponding eigenvalue governed by:

$$\mathbf{V}\_{\mathbf{n}}^{\prime\prime} + \left(\boldsymbol{\beta}\_{\mathbf{n}}^{2} - \frac{\mathbf{u}\_{0}^{\prime}}{\mathbf{u}\_{0}}\right) \mathbf{V}\_{\mathbf{n}} = \mathbf{0}.\tag{59}$$

$$\mathbf{V\_n(0)} = \mathbf{0},\tag{60}$$

$$\mathbf{V\_n(1)} = \mathbf{0}.\tag{61}$$

The solution of (59)–(61) is obtained numerically subject to the additional normalization condition V 0 n (0) = 1. Figure 6 depicts the profiles of the shape functions for the first five modes for Q = 0.35. We observe that β<sup>n</sup> increases with n (see Appendix B). For large n, problem (59)–(61) tends to a Sturm-Liouville problem with eigenvalues β<sup>n</sup> ∼ nπ and trigonometric shape functions Vn(z) ∼ sin(nπz) nπ . The profiles in Figure 6 clearly acquire the trigonometric character for larger n values. The profiles for other Q values present the same trend.

*Fluids* **2020**, *5*, x 25 of 51

**Figure 6.** Eigenfunction distributions for the first five modes at Q = 0.35. **Figure 6.** Eigenfunction distributions for the first five modes at Q = 0.35.

The stream function, velocity components and pressure to O Re−<sup>1</sup> are then given by

$$\Psi(\mathbf{x} < \mathbf{0}, \mathbf{z}) = (\mathbf{1} - \mathbf{G})\frac{\mathbf{z}^2}{2} + \mathbf{G}\frac{\mathbf{z}^3}{3} + \mathbf{Re}^{-1}\sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{e}^{\beta\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}(\mathbf{z}),\tag{62}$$

$$\mathbf{u}(\mathbf{x} < \mathbf{0}, \mathbf{z}) = (1 - \mathbf{G})\mathbf{z} + \mathbf{G}\mathbf{z}^2 + \mathbf{Re}^{-1} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\beta\_{\mathbf{n}}} \mathbf{e}^{\beta\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}(\mathbf{z}),\tag{63}$$

$$\mathbf{w}(\mathbf{x} < \mathbf{0}, \mathbf{z}) = -\mathrm{Re}^{-1} \sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{e}^{\beta\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}(\mathbf{z}) \,\tag{64}$$

n 1 n

 

and

<sup>A</sup> u x 0, z 0 Re e 0

$$\mathbf{p}(\mathbf{x}\leq \mathbf{0}, \mathbf{z}) = \mathbf{R}\mathbf{e}^{-1} \left[ \mathbf{2Gx} - \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{e}^{\mathfrak{P}\_{\mathbf{n}}\mathbf{x}} \left( \left[ (\mathbf{1} - \mathbf{G})\mathbf{z} + \mathbf{G}\mathbf{z}^{2} \right] \mathbf{V}'\_{\mathbf{n}} - \left[ (\mathbf{1} - \mathbf{G}) + \mathbf{2G}\mathbf{z} \right] \mathbf{V}\_{\mathbf{n}} \right) \right]. \tag{65}$$

n 1 1 2 <sup>n</sup> <sup>x</sup> <sup>n</sup> n n n 1 n <sup>A</sup> p x 0, z Re 2Gx e 1 G z Gz V 1 G 2Gz V . (65) We observe that the pressure at the stationary die or blade (z = 0) is simply Rep x 0,z 0 2Gx , corresponding to fully-developed flow. Additional correction is determined We observe that the pressure at the stationary die or blade (z = 0) is simply Rep(x ≤ 0, z = 0) = 2Gx, corresponding to fully-developed flow. Additional correction is determined when we examine the boundary layer near the stationary wall in Section 5. In this regard, we emphasize that the solution in the core layer is predominantly inviscid rotational. Consequently, the flow field given by (62)–(65) is not expected to satisfy adherence at the moving substrate (z = 1) and the stationary slot or blade wall (z = 0). This becomes particularly evident when examining the velocity at the walls: <sup>u</sup>(<sup>x</sup> <sup>&</sup>lt; 0, z <sup>=</sup> <sup>0</sup>) <sup>=</sup> Re−<sup>1</sup> <sup>P</sup><sup>∞</sup> n=1 An β<sup>n</sup> e <sup>β</sup>n<sup>x</sup> , <sup>0</sup> and <sup>u</sup>(<sup>x</sup> <sup>&</sup>lt; 0, z <sup>=</sup> <sup>1</sup>) <sup>=</sup> <sup>1</sup> <sup>+</sup> Re−<sup>1</sup> <sup>P</sup><sup>∞</sup> n=1 An β<sup>n</sup> e <sup>β</sup>nxV 0 n (z = 1) , 1. This signals the presence of the slip layers III and IV in Figure 1, which will be dealt with in Sections 5 and 6, respectively.

when we examine the boundary layer near the stationary wall in Section 5. In this regard, we emphasize that the solution in the core layer is predominantly inviscid rotational. Consequently, the flow field given by (62)–(65) is not expected to satisfy adherence at the moving substrate (z = 1) and the stationary slot or blade wall (z = 0). This becomes particularly evident when examining the The coefficients A<sup>n</sup> are determined by matching the inside and outside flows at the exit, x = 0, and will therefore be determined after the outside flow is examined next. Thus, for x > 0, the solution of problem (54)–(57) takes the following form:

velocity at the walls: <sup>1</sup> <sup>n</sup> <sup>x</sup> <sup>n</sup>

*Fluids* **2020**, *5*, 180 free-surface layer. It satisfies the following equation and boundary conditions:

$$\overline{\mathbf{w}}(\mathbf{x}>0,\mathbf{z}) = 2\frac{\mathbf{G}}{1-\mathbf{G}}\mathbf{V}\_0(\mathbf{z}) + \sum\_{n=1}^{\infty} \mathbf{A}\_n \mathbf{e}^{-\beta\_n n} \mathbf{V}\_n(\mathbf{z}).\tag{66}$$

The additional term V0(z) contributes to the matching of the core flow with the flow in the free-surface layer. It satisfies the following equation and boundary conditions: <sup>0</sup> z 1 2k V z 2kz 1 kz ln 1 2kz z 1 kz 2k ln 1 k 

*Fluids* **2020**, *5*, x 26 of 51

. This signals the presence of the slip layers III and IV in

The coefficients An are determined by matching the inside and outside flows at the exit, x = 0,

and will therefore be determined after the outside flow is examined next. Thus, for x > 0, the solution

<sup>x</sup> <sup>n</sup> 0 nn

The additional term V z <sup>0</sup> contributes to the matching of the core flow with the flow in the

<sup>G</sup> w x 0, z 2 V z A e V z 1 G

n 1

. (66)

$$\mathbf{V}\_0'' - \frac{\mathbf{u}\_0''}{\mathbf{u}\_0} \mathbf{V}\_0 = \mathbf{0}, \ \mathbf{V}\_0(0) = 1, \ \mathbf{V}\_0(1) = 0,\tag{67}$$

which admits an analytical solution: where <sup>G</sup> <sup>k</sup> . It is useful to examine the first contribution in (66) as it is the term that survives with

1 G

<sup>1</sup> <sup>n</sup> <sup>x</sup> <sup>n</sup> <sup>n</sup>

n 1 n <sup>A</sup> u x 0, z 1 1 Re e V z 1 1 

Figure 1, which will be dealt with in Sections 5 and 6, respectively.

of problem (54)–(57) takes the following form:

$$\mathbf{V}\_{0}(\mathbf{z}) = 2\mathbf{k}\mathbf{z}(1+\mathbf{k}\mathbf{z})\ln\left|\frac{\mathbf{z}}{1+\mathbf{k}\mathbf{z}}\right| + 1 + 2\mathbf{k}\mathbf{z} + \mathbf{z}(1+\mathbf{k}\mathbf{z})\left[2\mathbf{k}\ln|1+\mathbf{k}| - \frac{1+2\mathbf{k}}{1+\mathbf{k}}\right] \tag{68}$$

where k = <sup>G</sup> 1−G . It is useful to examine the first contribution in (66) as it is the term that survives with increasing downstream distance. The profiles of <sup>G</sup> <sup>1</sup>−<sup>G</sup> <sup>V</sup>0(z) are shown in Figure <sup>7</sup> for 0.4 <sup>≤</sup> <sup>Q</sup> <sup>≤</sup> 0.8. A couple of interesting distinctions can be made between slot coating (Q ≤ 0.5) and blade coating (<sup>Q</sup> <sup>≥</sup> 0.5). For slot coating, <sup>G</sup> <sup>1</sup>−<sup>G</sup> <sup>V</sup>0(z) decreases monotonically with z and remains positive. In contrast, it exhibits a minimum and remains negative for blade coating. This behavior suggests that the vertical flow is predominantly towards the moving substrate in slot coating, and away from it in blade coating. A couple of interesting distinctions can be made between slot coating Q 0.5 and blade coating Q 0.5 . For slot coating, <sup>0</sup> <sup>G</sup> V z 1 G decreases monotonically with z and remains positive. In contrast, it exhibits a minimum and remains negative for blade coating. This behavior suggests that the vertical flow is predominantly towards the moving substrate in slot coating, and away from it in blade coating.

**Figure 7.** Influence of flow rate on the profile of the dominant term in (66) or V0(z).

As mentioned above, the coefficients A<sup>n</sup> are found by matching the flow at the exit. For this purpose, equating the expressions (58) and (66) at x = 0 yields

$$\text{V}\_0(\mathbf{z}) = -\frac{1-\mathbf{G}}{\mathbf{G}} \sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_\mathbf{n} \mathbf{V}\_\mathbf{n}(\mathbf{z}) \,\tag{69}$$

which is a spectral representation of V0(z) in terms of the orthogonal shape functions Vn(z). Upon projecting and noting that <sup>h</sup>VnV0<sup>i</sup> <sup>=</sup> <sup>1</sup> β2 n , where the brackets < > denote integration over the interval 0 ≤ z ≤ 1, we have

$$\mathbf{A\_{n}} = -\left(\frac{\mathbf{G}}{1-\mathbf{G}}\right)\frac{\langle \mathbf{V\_{n}}\mathbf{V\_{0}}\rangle}{\left\langle \mathbf{V\_{n}^{2}}\right\rangle} = -\left(\frac{\mathbf{G}}{1-\mathbf{G}}\right)\frac{1}{\beta\_{\mathrm{n}}^{2}\langle \mathbf{v\_{n}^{2}}\rangle}.\tag{70}$$

A useful relation is obtained by first noting from (68) that V 0 0 (1) = −(1 − G). Consequently, upon differentiating (69) and evaluating it at z = 1, we have

$$\sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{V}\_{\mathbf{n}}'(1) = \mathbf{G}.\tag{71}$$

Finally, the flow field downstream of the die or blade exit become

$$\Psi(\mathbf{x} > \mathbf{0}, \mathbf{z}) = (\mathbf{1} - \mathbf{G})\mathbf{z}^2 + \mathbf{G}\frac{\mathbf{z}^3}{3} + \text{Re}^{-1} \left[ -2\frac{\mathbf{G}}{\mathbf{1} - \mathbf{G}}\mathbf{x}\mathbf{V}\_0(\mathbf{z}) + \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_\mathbf{n}}{\beta\_\mathbf{n}} \mathbf{e}^{-\beta\_\mathbf{n}\mathbf{x}}\mathbf{V}\_\mathbf{n}(\mathbf{z}) \right] \tag{72}$$

$$\mathbf{u}(\mathbf{x}>0,\mathbf{z}) = (1-\mathbf{G})\mathbf{z} + \mathbf{G}\mathbf{z}^2 + \mathbf{Re}^{-1} \left[ -2\frac{\mathbf{G}}{1-\mathbf{G}}\mathbf{x}\mathbf{V}\_0'(\mathbf{z}) + \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_\mathbf{n}}{\beta\_\mathbf{n}} \mathbf{e}^{-\beta\_\mathbf{n}\mathbf{x}}\mathbf{V}\_\mathbf{n}'(\mathbf{z}) \right].\tag{73}$$

$$\mathrm{cov}(\mathbf{x} > \mathbf{0}, \mathbf{z}) = \mathrm{Re}^{-1} \left[ 2 \frac{\mathbf{G}}{1 - \mathbf{G}} \mathbf{V}\_0(\mathbf{z}) + \sum\_{\mathbf{n} = 1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{e}^{-\beta\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}(\mathbf{z}) \right] \tag{74}$$

$$\mathbf{p}(\mathbf{x}\geq \mathbf{0}, \mathbf{z}) = -\mathrm{Re}^{-1} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathrm{A}\_{\mathrm{n}}}{\beta\_{\mathrm{n}}} \mathrm{e}^{-\beta\_{\mathrm{n}}\mathbf{x}} \left( \left[ (\mathbf{1} - \mathbf{G})\mathbf{z} + \mathbf{G}\mathbf{z}^{2} \right] \mathbf{V}'\_{\mathrm{n}} - \left[ (\mathbf{1} - \mathbf{G}) + 2\mathbf{G}\mathbf{z} \right] \mathbf{V}\_{\mathrm{n}} \right). \tag{75}$$

Here again, we see that the core solution does not ensure adherence at the moving substrate. This is clearly reflected in the value of the streamwise velocity at the moving substrate: u(x > 0, z = 1) = 1 + Re−<sup>1</sup> " 2Gx + P∞ n=1 An β<sup>n</sup> e <sup>−</sup>βnxV 0 n (1) # , 1, indicating the existence of a slip layer near the substrate, which will be discussed in Section 6. Interestingly, we can further see that the streamwise elongation rate component is also not zero since ux(x > 0, z = 1) = Re−<sup>1</sup> " 2G − P∞ n=1 Ane <sup>−</sup>βnxV 0 n (1) # , 0. In particular, using (71), we see that the exit ux(x = 0, z = 1) = GRe−<sup>1</sup> . These expressions will help estimate the thickness of the slip layer IV (see Section 6).

Unlike the velocity, expressions (65) and (75) for the core pressure turn out to be valid everywhere: in the core layer as well as at the stationary and the moving walls and near the free surface. The uniform validity of the pressure can be shown by determining the composite solution across the film [31,37]. Physically, the reason for the validity of the core pressure in the boundary layers is the result of the hydrostatic nature of the pressure in each layer and its invariability across each layer. Interestingly, while the core pressure given by (65) and (75) yields a correction to the Poiseuille level at any location, this correction is not felt at the stationary slot or blade wall to the current order. In fact, setting z = 0 in (65) and recalling that <sup>V</sup>n(0) <sup>=</sup> <sup>0</sup> from (60), we obtain the linear behavior: <sup>p</sup>(<sup>x</sup> <sup>≤</sup> 0, z <sup>=</sup> <sup>0</sup>) <sup>=</sup> 2Re−1Gx. The departure from the Poiseuille level will be established once we examine the pressure in the slip layer III near the stationary wall in the next Section 5. In contrast, the departure becomes increasingly evident for z > 0 as one approaches the moving substrate for a given x position. Indeed, evaluating the pressure at the moving substrate (75) gives

$$\text{Rep}(\mathbf{x} \le \mathbf{0}, \mathbf{z} = 1) = \mathbf{2Gx} - \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A\_{\mathbf{n}}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{e^{\beta\_{\mathbf{n}} \mathbf{x}} V\_{\mathbf{n}}'(1), \tag{76}$$

$$\text{Rep}(\mathbf{x}\geq 0, \mathbf{z}=1) = -\sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{F}\_{\mathbf{n}}} \mathbf{e}^{-\beta\_{\mathbf{n}}\mathbf{x}} \mathbf{V}\_{\mathbf{n}}'(1). \tag{77}$$

These expressions clearly signal the departure from the Poiseuille level. Figure 8 displays the pressure distribution along the moving substrate, typically illustrating the case for slot (Figure 8a) and blade (Figure 8b) coating. For each flow rate, the figure shows the linear behavior of the pressure starting far upstream, dictated by the 2Gx term in (76). Upon approaching the exit, the pressure experiences an exponential deviation, which is sustained beyond the exit as the pressure continues to decay and relaxes eventually to zero. The rate of relaxation is slower for a thinner film in slot coating and a thicker film in blade coating. A similar behavior is also predicted for the pressure distribution along the centerline of a symmetric planar jet as shown by Khayat [37] for the Newtonian limit (n = 1). Incidentally, the pressure continues to decrease for a shear-thinning fluid.

The departure from the Poiseuille level may be estimated by examining the pressure at the exit from (76)–(77), namely Rep(x = 0, z = 1) = − P∞ n=1 An β<sup>n</sup> V 0 n (1). This expression, however, is not particularly illuminating as the dependence on the flow rate is only implicit through the value of A<sup>n</sup> and β<sup>n</sup> , as well as V 0 n (1), as indicated in the table in Appendix B. In contrast, the influence of Q on the departure in the pressure gradient at the exit can be determined explicitly. Indeed, since the pressure gradient must match at the exit, then upon equating p<sup>x</sup> (x = 0, z = 1) from (76) and (77), we find that

$$\text{Rep}\_{\mathbf{x}}(\mathbf{x}=\mathbf{0}, \mathbf{z}=\mathbf{1}) = \sum\_{\mathbf{n}=1}^{\infty} \text{A}\_{\mathbf{n}} \text{V}\_{\mathbf{n}}'(1) = \mathbf{G} = \mathbf{3}(1 - \mathbf{2Q}), \tag{78}$$

which explicitly reflects the dependence of the departure from the constant upstream slope 2G = 6(1 − 2Q). Thus, the pressure gradient along the moving substrate evolves from 2G = 6(1 − 2Q) far upstream and reduces to G = 3(1 − 2Q) at the exit, which is also reflected Figure 8. Only the drag flow does not experience any change in the (zero) magnitude of the pressure gradient. Incidentally, (78) shows that (71) is equivalent to matching the pressure gradient at x = 0.

Interestingly, while the core pressure near the moving substrate exhibits a departure from the Poiseuille level, it retains the same behavior near the stationary die or blade. Another interesting observation worth making is the absence of transverse pressure gradient near both the stationary and moving walls. This is easily confirmed by showing from (65) and (75) that p<sup>z</sup> (x, z) vanishes at z = 0 an z = 1. The absence of a transverse pressure gradient suggests, in turn, the existence of a boundary or slip layer near each wall. The flow structure in each slip layer is examined in the next two sections.

#### **5. The Flow Near the Stationary Slot and Blade Walls (Slip Layer III)**

The flow structure in the boundary or slip layer III near the die or blade wall, the lower-wall layer, is examined in this section. In this layer, similar to the free-surface layer I, the transverse coordinate near the boundary will be taken as y = z = εη, where ε is the same small parameter used before. To examine the structure of layer III upstream of the channel exit, the near-wall coordinates are introduced as x = ξ and z = εη. Similar to the free-surface layer, matching with the core flow upstream of the channel exit shows that ψ ∼ (1 − G) y 2 <sup>2</sup> = (1 − G) ε 2η 2 2 as η → ∞. Therefore, ψ must be of O ε 2 close to the lower wall. In this case, the transformed Equation (4), along with conditions (10) and (11) become (dropping the overbar):

$$\begin{split} \psi\_{\eta} \Psi\_{\xi\eta} - \psi\_{\xi} \psi\_{\eta\eta} &= -\varepsilon^{2} \mathbf{p}\_{\xi} + \varepsilon^{2} \psi\_{\eta\eta\eta} + \varepsilon^{4} \psi\_{\xi,\xi\eta}, \\ \psi\_{\eta} \Psi\_{\xi\xi} - \psi\_{\xi} \psi\_{\xi\eta} &= \mathbf{p}\_{\eta} + \varepsilon^{2} \psi\_{\xi\eta\eta} + \varepsilon^{4} \psi\_{\xi\xi,\xi\zeta}, \\ \psi(\xi,\eta = 0) &= \psi\_{\eta}(\xi,\eta = 0) = 0, \\ \psi(\xi \to -\infty,\eta) &\sim (1 - \mathbf{G})\frac{\varepsilon^{2}\eta^{2}}{2} + \mathbf{G}\frac{\varepsilon^{3}\eta^{3}}{3}. \end{split} \tag{79}$$

flow does not experience any change in the (zero) magnitude of the pressure gradient. Incidentally,

Interestingly, while the core pressure near the moving substrate exhibits a departure from the Poiseuille level, it retains the same behavior near the stationary die or blade. Another interesting observation worth making is the absence of transverse pressure gradient near both the stationary and moving walls. This is easily confirmed by showing from (65) and (75) that <sup>z</sup> p x,z vanishes at z =

(78) shows that (71) is equivalent to matching the pressure gradient at x = 0.

**Figure 8.** Influence of the flow rate on the pressure distribution along the moving substrate for (**a**) Q 0.5 and (**b**) Q 0.5 . Drag flow corresponds to Q = 0.5. **Figure 8.** Influence of the flow rate on the pressure distribution along the moving substrate for (**a**) Q ≤ 0.5 and (**b**) Q ≥ 0.5. Drag flow corresponds to Q = 0.5.

The streamwise and transverse velocity components become u = <sup>1</sup> <sup>ε</sup>ψ<sup>η</sup> and w = −ψξ, respectively. Noting from the expressions (79) that p is of order ε 2 , the flow field expansion reads

$$
\psi(\xi,\eta) = \varepsilon^2 \Psi\_2(\xi,\eta) + \varepsilon^3 \Psi\_3(\xi,\eta) + \cdots,\tag{80}
$$

$$\mathbf{p}(\boldsymbol{\xi},\boldsymbol{\eta}) = \varepsilon^2 \mathbf{P}\_2(\boldsymbol{\xi},\boldsymbol{\eta}) + \varepsilon^3 \mathbf{P}\_3(\boldsymbol{\xi},\boldsymbol{\eta}) + \cdots \ . \tag{81}$$

On inserting (80) and (81) in (79), we obtain a hierarchy of problems to different orders, each problem requiring only one boundary condition in the x direction. The problem is therefore well posed by imposing the far-upstream flow as condition. Consequently, the matching of the upstream flow with the flow in the free-surface layer II at the channel exit (x = 0) is not necessary. Moreover, the rescaling in the streamwise direction is not required unless one wants to capture the flow structure very close to the exit, and the distance x is assumed to remain at least of O(1). In contrast to the flow in a constricted (dilated) channel [51], where the streamwise direction x is rescaled in terms of the inverse power of the indentation slope (and the Reynolds number), the flow is not captured very close to the origin in the present formulation.

We determine the additional boundary conditions to solve the problem (79) by matching between the lower-wall and the core layers. Proceeding as in Section 2 and Appendix A, we find that Ψ<sup>2</sup> ∼ (1 − G) η 2 2 and Ψ<sup>3</sup> ∼ G η 3 3 for large η. It is not difficult to show that the solution is trivial in both cases. Thus,

$$\Psi\_2(\xi,\eta) = (1-\mathcal{G})\frac{\eta^2}{2},\ P\_2(\xi,\eta) = 0,\tag{82}$$

$$\mathbb{V}\_{\mathfrak{P}}(\xi,\mathfrak{v}) = \mathrm{G}\frac{\mathfrak{v}^{3}}{\mathfrak{Z}}\mathsf{P}\_{\mathfrak{Z}}(\xi,\mathfrak{v}) = \mathfrak{Z}\mathsf{G}\xi.\tag{83}$$

Thus, up to O ε <sup>3</sup> = Re−<sup>1</sup> , the flow retains its Couette–Poiseuille character near the stationary die or blade wall. The correction is established when considering the next order. The equations for Ψ<sup>4</sup> and P<sup>4</sup> are

$$
\Psi\_{2\eta} \Psi\_{4\xi,\eta} - \Psi\_{4\xi} \Psi\_{2\eta\eta} = -\mathcal{P}\_{4\xi} + \Psi\_{4\eta\eta\eta\prime} \mathcal{P}\_{4\eta} = 0. \tag{84}
$$

Matching with the core flow by equating E3H4ψ = H4E3ψ yields the desired matching condition Ψ4(ξ, η → ∞) ∼ ηψ3z(x, z = 0) where, upon recalling that V 0 n (z = 0) = 1, we obtain ψ3z(ξ, z = 0) = P∞ n=1 An β<sup>n</sup> e <sup>β</sup>n<sup>ξ</sup> from (62). By using (82) and eliminating the pressure from (84), the problem for Ψ4(ξ, η) becomes,

$$\begin{aligned} (\mathbf{1} - \mathbf{G})\boldsymbol{\eta}\Psi\_{4\xi\boldsymbol{\eta}\eta\boldsymbol{\eta}} &= \Psi\_{4\eta\boldsymbol{\eta}\eta\boldsymbol{\eta}\boldsymbol{\eta}}\\ \Psi\_{4}(\boldsymbol{\xi},\boldsymbol{\eta}=\mathbf{0}) = \Psi\_{4\eta}(\boldsymbol{\xi},\boldsymbol{\eta}=\mathbf{0}) &= \mathbf{0}, \qquad \Psi\_{4}(\boldsymbol{\xi}\rightarrow-\boldsymbol{\infty},\boldsymbol{\eta}) \sim \mathbf{0},\\ \Psi\_{4}(\boldsymbol{\xi},\boldsymbol{\eta}\rightarrow\boldsymbol{\infty}) &\sim \mathbf{\eta}\sum\_{\mathbf{n}=1}^{\infty}\frac{\mathbf{A}\_{\mathbf{n}}}{\mathbf{B}\_{\mathbf{n}}}\mathbf{e}^{\beta\_{\mathbf{n}}\cdot\boldsymbol{\xi}}. \end{aligned} \tag{85}$$

The solution of problem (85) is not trivial, signaling the departure from the C-P flow. We seek a similarity solution of the form Ψ4(ξ, η) = P∞ n=1 An β<sup>n</sup> e <sup>β</sup>nξFn(η), where the functions Fn(η) are governed by

$$\mathbf{F\_{n}^{\prime\prime}} = (1 - \mathbf{G})\boldsymbol{\eta}\boldsymbol{\beta\_{n}}\mathbf{F\_{n}^{\prime}}$$

$$\mathbf{F\_{n}(0)} = \mathbf{F\_{n}^{\prime}(0)} = \mathbf{0}, \qquad \mathbf{F\_{n}^{\prime}(\infty)} \sim \mathbf{1}, \qquad \mathbf{F\_{n}^{\prime}(\infty)} \sim \mathbf{0}.\tag{86}$$

Although an analytical solution is possible in terms of Airy functions, problem (86) can be more conveniently solved numerically. The problem is further simplified by eliminating the explicit dependence on G and β<sup>n</sup> by introducing the transformation η → (1 − G) 1/3 β 1/3 <sup>n</sup> η and F<sup>n</sup> → (1 − G) 1/3 β 1/3 <sup>n</sup> F<sup>n</sup> .

Figure 9 illustrates the profiles of F<sup>n</sup> and its derivatives plotted against η, reflecting a monotonic behavior of Fn(η) and its derivatives. Of particular physical significance are the values F 00 n (η = 0) = c3(1 − G) 1/3 β 1/3 n and F<sup>000</sup> n (η = 0) = c4(1 − G) 2/3 β 2/3 n , which, as we shall see, are directly related to the shear stress at the wall and pressure, respectively. Here c<sup>3</sup> = 1.0651 and c<sup>4</sup> = −0.7765.

*Fluids* **2020**, *5*, x 32 of 51

**Figure 9.** The function Fn and its derivatives plotted against . Shown are: 1/3 1/3 n n 1G F as **Figure 9.** The function F<sup>n</sup> and its derivatives plotted against η. Shown are: (1 − G) 1/3β 1/3 <sup>n</sup> F<sup>n</sup> as curve (1), F0 <sup>n</sup> as curve (2), (1 − G) <sup>−</sup>1/3β −1/3 <sup>n</sup> F 00 n as curve (3) and (1 − G) <sup>−</sup>2/3β −2/3 <sup>n</sup> F 000 n as curve (4).

curve (1), Fn as curve (2), 1/3 1/3 n n 1G F as curve (3) and 2/3 2/3 n n 1G F as curve (4). To the current order, we determine the flow field in the slip layer III by inserting (82) and (83) as well as Ψ<sup>4</sup> into (80), yielding

$$\psi(\mathbf{x}, \mathbf{z}) = \frac{1 - \mathbf{G}}{2} \mathbf{z}^2 + \frac{\mathbf{G}}{3} \mathbf{z}^3 + \text{Re}^{-4/3} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\beta\_{\mathbf{n}}} \mathbf{e}^{\beta\_{\mathbf{n}} \times} \mathbf{F}\_{\mathbf{n}} \{ \text{Re}^{1/3} \mathbf{z} \}. \tag{87}$$

 2 3 4/3 1/3 <sup>n</sup> <sup>x</sup> <sup>n</sup> 1G G <sup>A</sup> x,z z z Re e F Re z 2 3 . (87) Of particular interest are the pressure and the shear stress at the wall:

$$\text{Rep}(\mathbf{x}) = 2\mathbf{G}\mathbf{x} + \text{Re}^{-1/3}\mathbf{c}\_4 (1 - \mathbf{G})^{2/3} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_\mathbf{n}}{\mathfrak{F}\_\mathbf{n}^{1/3}} \mathbf{e}^{\mathfrak{G}\_\mathbf{n}\mathbf{x}},\tag{88}$$

n

$$\tau\_{\infty}(\mathbf{x}, \mathbf{z} = 0) = 1 - \mathbf{G} + \mathbf{c}\_{3} (1 - \mathbf{G})^{1/3} \mathbf{R} \mathbf{e}^{-2/3} \sum\_{\mathbf{n} = 1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{F}\_{\mathbf{n}}^{2/3}} \mathbf{e}^{\mathfrak{B}\_{\mathbf{n}} \mathbf{x}}.\tag{89}$$

n 1 n 1/3 2/3 <sup>n</sup> xn xz 3 2/3 n 1 n <sup>A</sup> x,z 0 1 G c 1 G Re e . (89) We recall that the pressure is hydrostatic across the slip layer III, which is reflected by the absence of the z dependence in (88). Figure 10 illustrates the influence of the flow rate on the pressure (Figure 10a) and the wall shear stress (Figure 10b) at Re = 10. The range Q < 0.5 is taken to correspond typically to slot coating, and Q = 0.5 corresponds to drag flow. In this case, the adverse pressure gradient in Figure 10a opposes the action of the forward translation of the substrate, exhibiting the constant positive gradient 2G 6 1 2Q far upstream. However, the pressure rises rather rapidly as the flow approaches the die exit after becoming positive at a point that depends on the flow rate. However, this dependence is not monotonic. As Q departs (decreases) from the drag-flow level, Q = 0.5, where there is no change in the pressure, the pressure rises at the exit, but only to reach a maximum around Q = 0.4 and drops again as illustrated for Q = 0.35. This non-monotonicity is a consequence of flow separation, which will be discussed shortly. More generally, the rise in pressure We recall that the pressure is hydrostatic across the slip layer III, which is reflected by the absence of the z dependence in (88). Figure 10 illustrates the influence of the flow rate on the pressure (Figure 10a) and the wall shear stress (Figure 10b) at Re = 10. The range Q < 0.5 is taken to correspond typically to slot coating, and Q = 0.5 corresponds to drag flow. In this case, the adverse pressure gradient in Figure 10a opposes the action of the forward translation of the substrate, exhibiting the constant positive gradient 2G = 6(1 − 2Q) far upstream. However, the pressure rises rather rapidly as the flow approaches the die exit after becoming positive at a point that depends on the flow rate. However, this dependence is not monotonic. As Q departs (decreases) from the drag-flow level, Q = 0.5, where there is no change in the pressure, the pressure rises at the exit, but only to reach a maximum around Q = 0.4 and drops again as illustrated for Q = 0.35. This non-monotonicity is a consequence of flow separation, which will be discussed shortly. More generally, the rise in pressure depends also on inertia; its behavior is estimated from (88) to be Rep(<sup>x</sup> <sup>=</sup> 0, z) ∼ −Re−1/3 . We observe that the explicit dependence on Q is not evident from (88) since the coefficients A<sup>n</sup> and eigenvalues β<sup>n</sup> also depend on Q (or G as shown in Appendix B). The shear stress plots in Figure 10b show that the stress experiences a drop near the exit. We have also included the shear stress curves that correspond to the Couette–Poiseuille level for each flow rate as asymptotes, which roughly help locate the inception of the slip layer. The rise in adverse pressure above

the Poiseuille level (Figure 10a) is caused by the flow acceleration as the film approaches the exit; the flow converges sharply as the film contracts at the exit (see Figure 2). The rise in the pressure is accompanied by a drop in the shear stress, resulting in a loss in forward flow momentum; the increase in pressure in the direction of the flow is, of course, akin to an increase in the potential energy of the fluid, leading to a reduced kinetic energy and a deceleration of the fluid. 10a) is caused by the flow acceleration as the film approaches the exit; the flow converges sharply as the film contracts at the exit (see Figure 2). The rise in the pressure is accompanied by a drop in the shear stress, resulting in a loss in forward flow momentum; the increase in pressure in the direction of the flow is, of course, akin to an increase in the potential energy of the fluid, leading to a reduced kinetic energy and a deceleration of the fluid.

*Fluids* **2020**, *5*, x 33 of 51

that the explicit dependence on Q is not evident from (88) since the coefficients An and eigenvalues

n also depend on Q (or G as shown in Appendix B). The shear stress plots in Figure 10b show that the stress experiences a drop near the exit. We have also included the shear stress curves that

locate the inception of the slip layer. The rise in adverse pressure above the Poiseuille level (Figure

**Figure 10.** Pressure (**a**) and shear stress (**b**) distributions for slot coating along the stationary die wall for Re = 10 and various coating thicknesses (flow rates). Also added are the Couette–Poiseuille levels **Figure 10.** Pressure (**a**) and shear stress (**b**) distributions for slot coating along the stationary die wall for Re = 10 and various coating thicknesses (flow rates). Also added are the Couette–Poiseuille levels in (**b**) to show the leading edge of the slip layer III far upstream.

in (**b**) to show the leading edge of the slip layer III far upstream.

We observe that the flow in the wall layer is slower than in the core layer, and therefore expect a greater influence of the increasing pressure gradient. This seems to be particularly the case for a thinner coating film, here typically illustrated by Q = 0.35 for slot coating in Figure 10. The adverse pressure gradient is sufficiently large for the shear stress to vanish and a separation to eventually occur, with flow reversal occurring as it separates from the die wall.

We consider next the case of blade coating, which is typically illustrated in Figure 11, where we display the pressure (Figure 11a) and the wall shear stress (Figure 11b) for Q > 0.5 and Re = 10. The drag flow (Q = 0.5) is again included for reference. The pressure at the blade exit drops below atmospheric for the same reason as the rise in pressure in slot coating, namely as a result of the meniscus curvature and the film acceleration while moving along curved streamlines (see Figure 4a). For a higher flow rate, the fall in pressure is sharper, with a corresponding sharper rise in the wall shear stress (Figure 11b). Thus, a thicker film has to adjust more rapidly in height and velocity as Figure 2 indicates. Although the profiles in Figure 11 somewhat mirror those in slot coating, there are two important distinctions to observe. While the deviation from the Poiseuille level is non-monotonic with respect to Q for slot coating, the response is monotonic for blade coating (Figure 11a). This monotonicity is also reflected in the location of the inception of the slip layer, which coincides with the location of departure of the shear stress from the Couette-Poiseuille level, as illustrated in Figure 11b. More importantly, in contrast to slot coating, there is no possibility for separation in blade coating since there is no adverse pressure gradient (Figure 11a) and the wall shear stress remains positive (Figure 11b). Finally, the drop in pressure below atmospheric has been reported in the literature, which we will elucidate further next.

In their computational analysis of high-speed blade coating, Iliopoulos and Scriven [17] reported on pressure drop below atmospheric, similar to the drop predicted by the present formulation (Figure 11a). Direct full quantitative comparison is difficult since Iliopoulos and Scriven included the effect of forces such as shear thinning, the elastic deformation of the substrate and blade, as well as the wear on the blade from particle collisions, which are not accounted for in our study. Nevertheless, we will see that the comparison reveals important fundamental agreement and discrepancies when inertia is involved. In any case, the effect of the various additional forces do not seem to make any significant qualitative and quantitative difference, which can be confirmed by referring to figures of Iliopoulos and Scriven [17]. These figures indicate that the pressure upstream of the blade rises gradually from ambient to a sharp peak at the entrance to the blade region. The adverse pressure gradient causes the deflection of the excess liquid. The pressure then decreases linearly across the channel between the blade and the moving substrate, and drops below atmospheric at the blade exit. As mentioned earlier, this further drop is caused by the meniscus curvature and the film acceleration while moving along curved streamlines. The pressure then rises gradually back to atmospheric level as the film tends to uniform (plug) flow conditions.

*Fluids* **2020**, *5*, x 35 of 51

**Figure 11.** Pressure (**a**) and shear stress (**b**) distributions along the stationary blade or die wall for Re = 10 and various coating thicknesses Q 1/2 . Also added are the Couette–Poiseuille levels in (**b**) **Figure 11.** Pressure (**a**) and shear stress (**b**) distributions along the stationary blade or die wall for Re = 10 and various coating thicknesses (Q ≥ 1/2). Also added are the Couette–Poiseuille levels in (**b**) to show the leading edge of the wall layer far upstream.

to show the leading edge of the wall layer far upstream.

In their computational analysis of high-speed blade coating, Iliopoulos and Scriven [17] reported on pressure drop below atmospheric, similar to the drop predicted by the present formulation (Figure 11a). Direct full quantitative comparison is difficult since Iliopoulos and Scriven included the effect of forces such as shear thinning, the elastic deformation of the substrate and blade, as well as the wear on the blade from particle collisions, which are not accounted for in our study. Nevertheless, we will see that the comparison reveals important fundamental agreement and discrepancies when inertia is involved. In any case, the effect of the various additional forces do not seem to make any significant qualitative and quantitative difference, which can be confirmed by referring to figures of Iliopoulos and Scriven [17]. These figures indicate that the pressure upstream of the blade rises gradually from ambient to a sharp peak at the entrance to the blade region. The adverse pressure gradient causes the deflection of the excess liquid. The pressure then decreases linearly across the channel between the This behavior is also captured by the present theoretical analysis, which is depicted from the pressure profiles in Figure 12 for the three cases considered by Iliopoulos and Scriven [17], namely for Re = 1, 5 and 12 and corresponding flow rates of 0.58, 0.6 and 0.62, respectively. We have also included Figure 9 of Iliopoulos and Scriven as the top-right inset for reference, which depicts the pressure profiles for the trailing stiff blade geometry. We note that we used a much narrower streamwise range in Figure 12 than in the inset. In addition, our ambient pressure is 0 as opposed to 1 in the inset. Although the same flow rates are used by Iliopoulos and Scriven (inset), the present pressure gradients are higher. This is expected since the shear stress is higher for a Newtonian fluid. The flow rates used in Figure 12 were deduced from Figure 10 of Iliopoulos and Scriven [17] where the coating thickness is plotted against the Reynolds number. Incidentally, as they report: "Figure 10 suggests that for a stiff trailing blade with deformable substrate, raising the Reynolds number, or equivalently increasing the substrate speed, the liquid density, or decreasing the viscosity, leads to a thicker coating

film. For Re around 55, the coated thickness becomes equal to the gap between the blade and the substrate". This response is in line with the basic premise of the present theory which stipulates that the thickness increases with inertia, resulting from the drop in viscous effects responsible for the onset of the free-surface layer and the ensuing film contraction. In fact, the film thickness is equal to the channel gap as Re tends to infinity. In this case, the flow retains its CP profile across the exit region and further downstream. We emphasize that the CP is a solution of the Euler's equations. Another important agreement between the theoretical and numerical predictions is the dependence of the pressure drop below atmospheric on inertia. This is reflected in the inset at the bottom left of Figure 12, where we plot Rep(x = 0) against Re based on the three cases reported in the main figure. Clearly, the trend in the numerical and theoretical data points is essentially the same. The pressure at the exit increases rapidly in the small Re range and tapers (decaying to zero) for large Re. More precisely, expression (88) indicates that the pressure at the exit reduces to Rep(x = 0) = Re−1/3 c4(1 − G) 2/3 P<sup>∞</sup> n=1 An β 1/3 n , which suggests that the pressure drops below atmospheric is of order Re−1/3 . However, we note that if the flow rate (or G) varies with the Reynolds number as the data of Iliopoulos and Scriven [17] seem to suggest for a flexible substrate, then the dependence of the pressure drop on Re may differ slightly. We recall that A<sup>n</sup> and β<sup>n</sup> are G dependent (see Appendix B). In fact, it seems that Rep(0) behaves like Re−1/7 , a behavior which is surprisingly the same for both the theoretical and numerical predictions *Fluids* **2020**, *5*, x 37 of 51

**Figure 12.** Pressure profiles along the stationary blade (x < 0) and the free surface (x > 0) for Re = 1, 5 and 12 (Q = 0.58, 0.6 and 0.62). Inset on top right shows Figure 13 from Iliopoulos and Scriven [17], and inset on bottom left shows the values of Rep(0) against Re from current theory (squares) and the **Figure 12.** Pressure profiles along the stationary blade (x < 0) and the free surface (x > 0) for Re = 1, 5 and 12 (Q = 0.58, 0.6 and 0.62). Inset on top right shows Figure 13 from Iliopoulos and Scriven [17], and inset on bottom left shows the values of Rep(0) against Re from current theory (squares) and the numerical results of Iliopoulos and Scriven [17] (circles).

numerical results of Iliopoulos and Scriven [17] (circles).

**6. The Flow near the Moving Substrate (Slip Layer IV)** 

parameter in the problem, defined as Re or

<sup>3</sup> <sup>n</sup> xn <sup>n</sup>

p , ~ 2Gx e V 1

becomes (dropping the bar):

n 1 n A

proceeding as before, we see that in order to match the flow at the edge of the layer IV to the core layer I, the stream function and pressure must tend to , ~ and

, respectively, to lowest order in . In this case, for the inertial

and viscous terms to balance in the transformed momentum equations, we must have 1/2 . Therefore, 2 3 Re and 3/2 . In this case, the problem for , and p ,

To examine the structure of the slip layer IV near the moving substrate, we let y1z 0 . The scaling in the transverse direction is changed by writing y 0 , where is the small

1/ Re , where is to be determined.

*Fluids* **2020**, *5*, x 42 of 51

**Figure 13.** Influence of the flow rate on the shear stress distribution along the moving substrate at Re = 10 for (**a**) Q 1/2 and (**b**) Q 1/2 . Drag flow corresponds to Q 1/2 . Corresponding **Figure 13.** Influence of the flow rate on the shear stress distribution along the moving substrate at Re = 10 for (**a**) Q ≤ 1/2 and (**b**) Q ≥ 1/2. Drag flow corresponds to Q = 1/2. Corresponding fully-developed profiles are included in lighter lines.

fully-developed profiles are included in lighter lines. **7. Discussion and Concluding Remarks**  The planar laminar free-surface coating flow of a Newtonian fluid is investigated in the current study. The flow near and far from the channel exit is examined at moderate Reynolds and capillary numbers, subject to the substrate translation, and an adverse or a favorable constant pressure gradient applied far upstream of the channel exit as encountered in slot and blade coating flows. Figure 12 reveals, however, two fundamental features that do not seem to be captured by the numerical simulation (top-right inset). The first is the accelerated drop below atmospheric near the exit (main Figure 12) as opposed to the continuous linear drop in the inset. This means that the numerical simulation does not indicate the existence of any upstream influence, which is rather inaccurate given the elliptic nature of the governing equations. The second discrepancy concerns the stress singularity at the exit, which should inevitably be reflected in the pressure singularity as shown in figure but not

Although the flow far downstream is relatively simple to analyze as it becomes of the boundary-layer

in the inset. The pressure singularity is a consequence of elongational effect; near x = 0, the streamwise momentum equation reduces to the balance between the pressure gradient and the gradient of the excess normal stress, thus yielding Rep<sup>x</sup> = uxx. Consequently, the jump in u<sup>x</sup> from zero to a large positive value (see Figure 2b) leads to the jump in the pressure. The pressure singularity can manifest itself numerically in the form of spikes as reported in Figure 11b of Mitsoulis and Athanasopoulos [50]. We suspect that the singularity was smoothed over in the finite-element calculation of Iliopoulos and Scriven [17]. As we shall see next, the flow is very different in the slip layer IV along the moving substrate. We have already reported on the pressure profiles in Figure 8, which turn out to be smooth. The shear stress profiles require further development as we shall see in the next section.

#### **6. The Flow Near the Moving Substrate (Slip Layer IV)**

To examine the structure of the slip layer IV near the moving substrate, we let y = 1 − z > 0. The scaling in the transverse direction is changed by writing y = γη > 0, where γ is the small parameter in the problem, defined as γ = Re−<sup>β</sup> or Re = γ <sup>−</sup>1/β, where β is to be determined. Similar to the analysis of the free-surface layer II and the slip layer III, the following change of coordinates is introduced, namely x = ξ, z = 1 − γη. Letting ψ(ξ, η) = 1 − G <sup>6</sup> + ψ(ξ, η) and proceeding as before, we see that in order to match the flow at the edge of the layer IV to the core layer I, the stream function and pressure must tend to ψ(ξ, η → ∞) ∼ −γη and p(ξ, η → ∞) ∼ ε 3 " 2Gx − P∞ n=1 An β<sup>n</sup> e <sup>β</sup>nxV 0 n (1) # , respectively, to lowest order in γ. In this case, for the inertial and viscous terms to balance in the transformed momentum equations, we must have β = 1/2. Therefore, Re = γ <sup>−</sup><sup>2</sup> = ε <sup>−</sup><sup>3</sup> and γ = ε 3/2 . In this case, the problem for ψ(ξ, η) and p(ξ, η) becomes (dropping the bar):

$$\begin{split} \psi\_{\eta} \Psi\_{\xi,\eta} - \psi\_{\xi} \psi\_{\eta\eta} &= -\varepsilon^{3} \mathbf{p}\_{\xi} - \varepsilon^{3/2} \psi\_{\eta\eta\eta} - \varepsilon^{9/2} \psi\_{\xi,\xi,\eta}, \\ \psi\_{\eta} \Psi\_{\xi,\xi} - \psi\_{\xi} \psi\_{\xi,\eta} &= \mathbf{p}\_{\eta} - \varepsilon^{3/2} \psi\_{\xi,\eta\eta} - \varepsilon^{9/2} \psi\_{\xi,\xi,\xi,\eta} \\ \psi(\xi,\eta = 0) &= 0, \qquad \psi\_{\overline{\eta}}(\xi,\eta = 0) = -\varepsilon^{3/2}, \\ \psi(\xi \to -\infty, \eta) &\sim -\varepsilon^{3/2} \eta + \varepsilon^{3} \eta^{2} \Big(\frac{1+\mathcal{G}}{2}\Big) - \mathcal{G} \varepsilon^{9/2} \frac{\eta^{3}}{3}. \end{split} \tag{90}$$

Additional conditions for the stream function and the pressure are established from matching with the core flow. However, these conditions are not the same upstream (x < 0) and downstream (x > 0) of the exit. The two problems will be treated separately in each region. However, for any x, the expansions of the stream function and the pressure take the forms:

$$\Psi(\xi,\eta) = 1 - \frac{\mathcal{G}}{6} + \varepsilon^{3/2} \Psi\_2(\xi,\eta) + \varepsilon^3 \Psi\_3(\xi,\eta) + \varepsilon^{9/2} \Psi\_{9/2}(\xi,\eta) + \dots \tag{91}$$

$$\mathbf{p}(\xi,\mathfrak{n}) = \varepsilon^3 \mathbf{P}\_3(\xi,\mathfrak{n}) + \varepsilon^4 \mathbf{P}\_4(\xi,\mathfrak{n}) + \dots \tag{92}$$

Proceeding as in the previous sections, we find that the first two terms are contributions to the Couette-Poiseuille flow:

$$
\Psi\_{3/2}(\xi,\eta) = -\eta.\tag{93}
$$

$$\Psi\_3(\xi,\eta) = \left(\frac{1+G}{2}\right)\eta^2. \tag{94}$$

The correction occurs at the next order, O ε 6 , and the problem is governed by

$$\Psi\_{9/2\xi\eta} = \mathbf{P}\_{3\xi} + \Psi\_{9/2\eta\eta\eta} \tag{95}$$

$$P\_{\mathfrak{N}\mathfrak{h}} = 0,\tag{96}$$

$$\Psi\_{\mathcal{B}/2}(\xi\_{\prime}\eta = 0) = \Psi\_{\mathcal{B}/2\overline{\eta}}(\xi\_{\prime}\eta = 0) = 0,\tag{97}$$

$$\Psi\_{9/2}(\xi \to -\infty, \mathfrak{n}) \sim -\frac{\mathcal{G}}{3} \mathfrak{n}^3, \,\,\Psi\_{9/2}(\xi, \mathfrak{n} \to \infty) \sim -\frac{\mathcal{G}}{3} \mathfrak{n}^3. \tag{98}$$

We consider the solution of problem (95)–(98) separately upstream and downstream of the exit. For x < 0, the pressure remains dictated by its core value from (65):

$$\mathbf{p}\_3(\mathbf{x} \le \mathbf{0}, \boldsymbol{\eta} \to \infty) \sim 2\mathbf{G}\mathbf{x} - \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{e}^{\mathfrak{P}\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}'(1),\tag{99}$$

which, when inserted in (95), leads to the following equation for Ψ9/2 (ξ < 0, η):

$$\Psi\_{\mathfrak{P}/2\xi\eta} = -\sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{e}^{\mathfrak{P}\_{\mathbf{n}} \mathbf{x}} \mathbf{V}'\_{\mathbf{n}}(1) + \Psi\_{\mathfrak{P}/2\eta\eta\eta} + 2\mathbf{G}.\tag{100}$$

The solution may be written as Ψ9/2 (ξ < 0, η) = − G 3 η <sup>3</sup> + P∞ n=1 Ane <sup>β</sup>nxV 0 n (1)Gn(η), where the coefficients Gn(η) are governed by

$$\text{G}\_{\text{n}}^{\prime\prime\prime} - \beta\_{\text{n}} \text{G}\_{\text{n}}^{\prime} - 1 = 0, \text{ G}\_{\text{n}}(0) = \text{G}\_{\text{n}}^{\prime}(0), \text{ G}\_{\text{n}}^{\prime\prime\prime}(\infty) \sim 0. \tag{101}$$

This problem admits an analytical solution Gn(η) = −η + √ 1 β<sup>n</sup> 1 − e − √ βnη , and the stream function near the moving substrate becomes

$$\begin{split} \mathbb{Y}(\xi < 0, \eta) &= 1 - \frac{\mathbb{G}}{\mathbb{G}} - \mathrm{Re}^{-1/2}\eta + \mathrm{Re}^{-1} \Big( \frac{1+\mathcal{G}}{2} \Big) \eta^2 - \mathrm{Re}^{-3/2} \frac{\mathcal{G}}{\mathbb{G}} \eta^3 \\ &- \mathrm{Re}^{-3/2} \sum\_{n=1}^{\infty} \frac{\mathcal{A}\_n}{\overline{\beta}\_n} \mathrm{e}^{\beta\_n \chi} \mathcal{V}'\_n(1) \Big[ \eta - \frac{1}{\sqrt{\beta\_n}} \Big( 1 - \mathrm{e}^{-\sqrt{\beta\_n}\eta} \Big) \Big]. \end{split} \tag{102}$$

When expressed in terms of z, this expression leads to the following expression for the velocity and the shear stress near the moving substrate:

$$\mathbf{u}(\mathbf{x}<0,\mathbf{z}) = \mathbf{u}\_0(\mathbf{z}) + \mathrm{Re}^{-1} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{J}\_{\mathbf{n}}} \mathrm{e}^{\mathfrak{J}\_{\mathbf{n}} \times} \mathrm{V}\_{\mathbf{n}}'(1) \left(1 - \mathbf{e}^{-\sqrt{\mathfrak{J}\_{\mathbf{n}} \mathrm{Re}}(1-\mathbf{z})}\right) \tag{103}$$

$$\pi\_{\mathbf{x}\mathbf{x}}(\mathbf{x}<0,\mathbf{z}) = 1 - \mathbf{G} + 2\mathbf{G}\mathbf{z} - \mathrm{Re}^{-1/2} \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\sqrt{\mathfrak{B}\_{\mathbf{n}}}} \mathrm{e}^{\beta\_{\mathbf{n}}\mathbf{x}} V\_{\mathbf{n}}'(1) \mathrm{e}^{-\sqrt{\beta\_{\mathbf{n}}\mathbf{Re}(1-\mathbf{z})}}.\tag{104}$$

The exponential term that constitutes the main correction and ensures adherence at the moving substrate. Finally, the skin drag coefficient along the wall for x < 0 becomes <sup>τ</sup>xz(<sup>x</sup> <sup>&</sup>lt; 0, z <sup>=</sup> <sup>1</sup>) <sup>=</sup> <sup>1</sup> <sup>+</sup> <sup>G</sup> <sup>−</sup> Re−1/2 <sup>P</sup><sup>∞</sup> n=1 √ An β<sup>n</sup> e <sup>β</sup>nxV 0 n (1). This expression can be used to estimate the starting point of the slip layer IV, a point where the flow begins to deviate from the Couette-Poiseuille limit. Thus, the commencement of the slip layer would correspond, in practice, to a position x<sup>0</sup>

satisfying <sup>1</sup>+G−Re−1/2 <sup>P</sup><sup>∞</sup> n=1 √ An βn e <sup>β</sup>nx0V 0 n (1) 1+G equal to a small tolerance. Of course, the starting point depends

on Q and Re (see below). The slip layer IV should not be confused with the boundary layer examined by Carvalho and Kheshgi [9], which emanates at the dynamic contact line of the upstream meniscus in slot coating. Obviously, the upstream meniscus is not accounted for in the present formulation.

For x > 0, the problem is much more complicated, but remains mathematically manageable. The difficulty stems from the core pressure (65), which takes a different form than (99):

$$P\_3(\mathbf{x} > \mathbf{0}, \boldsymbol{\eta} \to \infty) \sim -\sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{e}^{-\mathfrak{P}\_{\mathbf{n}} \mathbf{x}} \mathbf{V}\_{\mathbf{n}}'(1),\tag{105}$$

which, when inserted in (95), leads to the following equation for Ψ9/2η(ξ > 0, η):

$$\Psi\_{9/2\xi,\eta} = \sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{e}^{-\beta\_{\mathbf{n}}\mathbf{x}} \mathbf{V}\_{\mathbf{n}}'(1) + \Psi\_{9/2\eta\eta\eta}.\tag{106}$$

In addition to conditions (97), Equation (106) must be solved subject to flow matching at the channel exit (x = 0):

$$\Psi\_{9/2\eta} \left( \xi = 0^+, \eta \right) = \Psi\_{9/2\eta} \left( \xi = 0^-, \eta \right) = -\mathbf{G}\eta^2 - \sum\_{\mathbf{n}=1}^{\infty} \frac{\mathbf{A}\_{\mathbf{n}}}{\mathfrak{P}\_{\mathbf{n}}} \mathbf{V}\_{\mathbf{n}}'(1) \left( 1 - \mathbf{e}^{-\sqrt{\beta\_{\mathbf{n}}}\eta} \right). \tag{107}$$

The solution may be conveniently written in the form

$$\Psi\_{9/2\overline{\eta}}(\xi\_{\nu}\overline{\eta}) = -\mathrm{G}\eta^2 - \sum\_{n=1}^{\infty} \mathrm{A}\_{\mathrm{n}} \mathrm{V}\_{\mathrm{n}}'(1) \left[ \mathrm{U}\_{\mathrm{n}}(\xi\_{\nu}\eta) + \frac{1}{\beta\_{\mathrm{n}}} \left( 1 - \mathrm{e}^{-\sqrt{\beta\_{\mathrm{n}}}\eta} \right) \right] \tag{108}$$

where the coefficients Un(ξ, η) are governed by the following problem

$$\begin{aligned} \mathbf{U\_{n\xi}} &= \mathbf{U\_{n\eta\eta}} - \mathbf{e^{-\sqrt{\beta\_n}\eta}} + \mathbf{2} - \mathbf{e^{-\beta\_n\chi}}, \\ \mathbf{U\_n(\xi, \eta = 0)} &= \mathbf{U\_n(\xi = 0^+, \eta)} = \mathbf{0}, \end{aligned} \tag{109}$$

which is of the transient heat conduction type with a source term, admitting the solution:

$$\begin{array}{lcl} \mathrm{U\_{n}}(\boldsymbol{\xi},\boldsymbol{\eta}) &= \int \limits\_{0}^{\xi} \frac{2-\mathrm{e}^{-\beta\boldsymbol{n}\_{s}}-\mathrm{e}^{-\sqrt{\beta\boldsymbol{n}\_{s}}}}{\sqrt{4\pi(\boldsymbol{\xi}-\boldsymbol{s})}} \Big[ \mathrm{exp} \Big( -\frac{\left(\boldsymbol{\eta}-\boldsymbol{y}\right)^{2}}{4\left(\boldsymbol{\xi}-\boldsymbol{s}\right)} \Big) - \mathrm{exp} \Big( -\frac{\left(\boldsymbol{\eta}+\boldsymbol{y}\right)^{2}}{4\left(\boldsymbol{\xi}-\boldsymbol{s}\right)} \Big) \Big] d\boldsymbol{y}d\boldsymbol{s} \\ &= \int \Big( 2-\mathrm{e}^{\beta\boldsymbol{n}\_{s}}(\boldsymbol{\xi}-\boldsymbol{t}) \Big) \mathrm{erf} \Big( -\frac{\eta}{\sqrt{4\boldsymbol{t}}} \Big( 1-\mathrm{e}^{\beta\boldsymbol{n}\_{s}\boldsymbol{t}} \Big) \mathrm{sinh} \Big( \sqrt{\beta\boldsymbol{n}\_{\boldsymbol{t}}} \eta \Big) \\ & - \frac{\boldsymbol{\xi}}{2} \int \mathrm{e}^{\beta\boldsymbol{n}\_{\boldsymbol{t}}} \Big[ \mathrm{e} \sqrt{\beta\boldsymbol{n}\_{\boldsymbol{t}}} \mathrm{erf} \Big( \frac{2\sqrt{\beta\boldsymbol{n}\_{\boldsymbol{t}}}\boldsymbol{t}+\eta}{2\sqrt{\boldsymbol{t}}} \Big) - \mathrm{e}^{-\sqrt{\beta\boldsymbol{n}\_{\boldsymbol{t}}}\eta} \mathrm{erf} \Big( \frac{2\sqrt{\beta\boldsymbol{n}\_{\boldsymbol{t}}}\boldsymbol{t}-\eta \Big) \Big] d\boldsymbol{t}. \end{array} \tag{110}$$

This, in turn, yields the following expression for the velocity to O Re−<sup>1</sup> near the moving substrate:

$$\mathbf{u}(\mathbf{x}>0,\mathbf{z}) = (1-\mathbf{G})\mathbf{z} + \mathbf{G}\mathbf{z}^2 + \mathbf{Re}^{-1} \sum\_{\mathbf{n}=1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{V}\_{\mathbf{n}}'(1) \left[ \frac{1}{\mathfrak{F}\_{\mathbf{n}}} \left( 1 - \mathbf{e}^{-\sqrt{\mathfrak{F}\_{\mathbf{n}}}\eta} \right) + \mathbf{U}\_{\mathbf{n}}(\xi,\eta) \right] \tag{111}$$

which is easily evaluated numerically. The shear stress or drag at the substrate takes the following form to O Re−1/2 :

$$\tau\_{\rm xx}(\mathbf{x} > 0, \mathbf{z} = 1) = 1 + \mathbf{G} - \mathbf{R}^{-1/2} \sum\_{\mathbf{n} = 1}^{\infty} \mathbf{A}\_{\mathbf{n}} \mathbf{V}\_{\mathbf{n}}'(1) \left[ 4 \sqrt{\frac{\mathbf{x}}{\pi}} + \frac{1 - \mathbf{e}^{\beta\_{\mathbf{n}} \cdot \mathbf{x}}}{\sqrt{\pi} \beta\_{\mathbf{n}}} + \frac{\mathbf{e}^{\beta\_{\mathbf{n}} \cdot \mathbf{x}} \text{erfc} \{ \sqrt{\beta\_{\mathbf{n}} \mathbf{x}} \}}{\sqrt{\beta\_{\mathbf{n}}}} \right]. \tag{112}$$

We recall again that the pressure is hydrostatic across the slip layer IV, and therefore remains as in Figure 8. Figure 13 illustrates the influence of the flow rate on the wall shear stress at Re = 10. The ranges Q < 0.5 and Q > 0.5 in Figure 13a,b are again taken to correspond typically to slot and blade coating, and Q = 0.5 corresponds to drag flow. We have also included the shear stress curves that correspond to the CP level for each flow rate as asymptotes, which roughly help locate the inception of the slip layer IV. For slot coating, Figure 13a indicates that the shear stress drops sharply from its CP level just before the exit, and can exhibit a minimum followed by a maximum before the sharp drop for a thin coating film as illustrated for Q = 0.35. For any Q < 0.5, Figure 13a suggests that the flow separates; the vanishing of the shear stress is accompanied by an adverse pressure (Figure 8a). It seems

that the separation always occurs at x < 0, at a location further upstream for the thinner coating film. In contrast, for blade coating, the shear stress plots in Figure 13b show that the stress experiences a rise near the exit. Although the stress can vanish for a thick coating film (Q > 0.6), there is no possibility of a separation given the favorable pressure gradient (see Figure 8b).

The rich dynamics and steep gradients exhibited by the shear stress for slot coating in Figure 13a, especially upstream of the exit, contrast sharply the smooth and mild behavior exhibited by the pressure in Figure 8a. This may seem at first physically unrealistic given the intimate coupling between stress and pressure gradients in the conservation equations. As a check, we examine the validity of the streamwise momentum equation along the wall for x < 0, which reduces to −Rep<sup>x</sup> + τxz,<sup>z</sup> = 0. Clearly, upon noting from (104) that τxz,<sup>z</sup> (x < 0, z = 1) = 2G − P∞ n=1 Ane <sup>β</sup>nxV 0 n (1), we deduce that Rep<sup>x</sup> (x < 0, z = 1) = 2G − P∞ n=1 Ane <sup>β</sup>nxV 0 n (1), which is the same result that is achieved upon evaluating the pressure gradient from (65) at z = 1. As to the transverse momentum equation, it reduces to −Rep<sup>z</sup> + τzx,<sup>x</sup> +τzz,<sup>z</sup> = 0 at z = 1. Noting that τzx,<sup>x</sup> +τzz,<sup>z</sup> = uzx + wzz = −wzz + wzz = 0 at the

moving substrate, this leaves p<sup>z</sup> = 0, which is the case to the current order as per (96). We thus confirm

that momentum is conserved in the streamwise and transverse directions.

## **7. Discussion and Concluding Remarks**

The planar laminar free-surface coating flow of a Newtonian fluid is investigated in the current study. The flow near and far from the channel exit is examined at moderate Reynolds and capillary numbers, subject to the substrate translation, and an adverse or a favorable constant pressure gradient applied far upstream of the channel exit as encountered in slot and blade coating flows. Although the flow far downstream is relatively simple to analyze as it becomes of the boundary-layer or the thin-film type, the treatment remains challenging in the vicinity of the exit. The method of matched asymptotic expansion is adopted to examine the influence of inertia and the applied pressure gradient on the shape of the free surface and the flow field. At the channel exit, a stress singularity occurs where the boundary condition changes from no-slip at the lower wall, to slip at the free surface, leading to the development of a boundary layer along the free surface. As a result, the flow domain consists of four different regions (Figure 1): the core layer I, the free-surface layer II, the slot or blade slip layer III, and the slip layer IV near the moving substrate. The layers II, III and IV are shear dominated, and the flow is obtained using a boundary-layer approach, but not all layers allow a similarity solution. In contrast, the core layer I is inviscid rotational where both shear and extensional flows are in balance. The small parameter in the problem turns out to be ε = Re−1/3 , based on the balance of inertia and viscous effects in the free-surface layer II, allowing the asymptotic development of the flow by expanding the flow field.

As the fluid emerges from the channel in the form of the jet, it experiences a drastic drop in the shear stress and rise in normal stress at the channel exit as it slips along the free surface. This type of singularity constitutes a major hurdle in a theoretical methodology, particularly in a computational approach. Moreover, there is a wealth of physical mechanisms and phenomena that cannot be easily captured by a numerical approach. The boundary-layer structures in layers II and III are investigated in detail. The similarity solutions obtained in the two layers do not require matching at the exit. Consequently, the presence of the singularity is circumvented, constituting a major advantage of the present formulation. The shape of the free surface is determined by matching the free-surface layer flow with the core flow outside the channel exit.

The inviscid character of the core layer I can be perplexing given the parabolic Couette-Poiseuille flow as the leading-order velocity profile. This is a pivotal point which embeds the major premise on which the present theory is built, and can be explained as follows. The Couette-Poiseuille flow is only recovered in the limit of infinite Reynolds number (ε → 0). In this limit, there is no viscous mechanism for the flow to change as it nears and traverses the exit. The fully developed Couette-Poiseuille velocity profile imposed far upstream retains its shape, and the free surface remains flat (horizontal). We therefore interpret the Couette-Poiseuille profile in the present context to correspond to inviscid

rotational flow. Mathematically, any fully-developed profile: u = u(z), w = 0, satisfies Euler's equations, which are the inviscid limit of the Navier–Stokes' equations for infinite Reynolds number. In the present problem, although the Couette-Poiseuille profile is derived by integrating the (viscous) Navier–Stokes' equations, this profile does indeed satisfy Euler's equations as well as all the boundary conditions at the solid walls and flat free surface in the inviscid limit. Physically, one expects that if the Reynolds number is large enough (excluding turbulence), one should observe a profile close to the parabolic profile, even downstream (such as the case of water out of a garden hose). This hypothesis, and the present approach, have been amply validated for an axisymmetric jet (see, for instance, Philippe and Dumargue [33]). Thus, in this inviscid limit, the parabolic (or any) profile imposed far upstream remains unchanged as the fluid emerges at the exit. There is simply no viscous mechanism for the shear stress to relax or (equivalently) for a boundary layer to form along the free surface. Consequently, the free surface remains horizontal. In fact, the shear stress vanishes everywhere in the inviscid limit (zero viscosity) regardless of the value of the velocity gradient. In this limit, the dynamic conditions (6) and (7) are identically satisfied.

The core layer of the slot or blade coating flow should, apparently, be similar to the channel flow with a fine constriction considered by Smith [51]. Both flows comprise a leading-order fully developed contribution and an inviscid rotational first-order correction. However, there is an important difference as a result of the presence of a free surface in coating flow. In coating, the fully developed flow prevails everywhere, including the region downstream of the die or blade exit, *only in the limit of infinite Reynolds* number. In this limit, the free surface remains flat. The situation is different in Smith's case. The Poiseuille flow is recovered in the limit of zero slope of the indentation *at any Reynolds number*, and therefore exists in that limit. In the absence of the constriction, the flow physically exists and is certainly viscous. Therefore, we tend to consider the Couette-Poiseuille flow in coating as inviscid as it is reached only in the limit of infinite Re, satisfying Euler's equations without any viscous mechanism capable of altering the Couette-Poiseuille profile as the flow traverses the exit. This is an unphysical limit flow. It is important to emphasize that whether the Couette-Poiseuille flow is labelled as inviscid or viscous has no consequence on the present development and results.

The free surface is found to always contract near the channel exit regardless of the level of inertia and direction of the applied pressure gradient. The film tends to slightly expand further downstream for a thick blade coating film (Figure 2). For Ca = O(1), we show that the shape of the meniscus is not affected by surface tension. However, surface tension appears to alter significantly the dynamics in the normal stress, especially for blade coating where it causes a maximum to emerge near the exit (Figure 3). Further downstream where the flow becomes of the boundary-layer type, we illustrate how the near-exit solution can be matched to a thin-film formulation (Figure 4). Alternatively, the near-exit solution can be used as initial condition for a computational approach, thus avoiding the incorporation of the singularity.

Given the inviscid character of the core layer I for any finite Reynolds number, the core solution does not satisfy adherence at the walls where slip layers III and IV emerge. The detailed flow structure is obtained in the slip layers. In particular, we find as a result of the vanishing wall shear stress and the presence of adverse pressure (Figure 10), that a separation may eventually occur for a thin film in slot coating. Separation appears to be even more unavoidable near the moving substrate for slot coating (Figure 13). Finally, the undershoot exhibited by the pressure at the blade is captured in Figure 12. Both the current and existing numerical results from the literature appear to suggest that the pressure drop below atmospheric grows like Re−1/7 . However, unlike existing numerical predictions that show a linear decay sustained all the way to the exit, the current predictions show a further nonlinear drop, thus reflecting an upstream influence not predicted when using a numerical approach.

We would like to conclude this section by discussing two important issues that are not addressed in the present study and are worth elaborating on for future consideration. The first being the case of small-flow rate (Q < 1/3) in slot coating. The present formulation cannot handle such a situation, which involves bending of the free surface as it invades the region between the channel walls upstream

of the exit. Consequently, the free-surface height experiences a significant departure from the infinite-Re limit of a flat surface, which violates the validity of expansion (19). We suspect that this would be an entirely different flow regime where the pressure gradient may no longer be assumed small relative to the Reynolds number. We also suspect that the flow in the immediate vicinity of the singularity with (approximate) solution to the full Navier–Stokes equations must be considered, which leads us to the second issue.

In fact, the second important issue concerns the flow very near the singularity. Much of twentieth-century boundary-layer theory was concerned with two problems, one the external flow over a plate (either infinite or semi-infinite length), the second flow in a slightly deformed channel by the presence of a mild constriction placed on the wall(s). The problem considered in this paper combines elements of both channel flow (through an upstream oncoming Couette–Poiseuille flow) with aspects of flow past a finite plate (through having the blade or die wall of the channel terminate with a free surface flow downstream of that point, so having elements of a wake-like flow). The present solution remains incomplete since the flow structure very near the blade or die exit is essentially left out. Only an approximation for the flow near the free surface downstream of the blade or die is derived. Since this boundary-layer solution (derived in Section 3) is only an approximate solution of the full Navier–Stokes equations, higher-order corrections remain envisageable.

The present development has a close parallel with the solution near the trailing point of a finite plate as proposed by Goldstein [47], who was able to develop approximate solutions to the wake, with the approximations including a singularity at the trailing edge. As in the present coating flow, the solution was developed in powers of x 1/3 , where x was the distance from the trailing edge, and the singularity was in effect an infinite transverse velocity (order x <sup>−</sup>2/3 as <sup>x</sup> <sup>→</sup> 0) and a singularity in the pressure and stress from x < 0 to x > 0. The singularity was circumvented, and Goldstein was able to analyze the leading order far wake flow. It took three decades before the asymptotic structure was properly explained, when Stewartson [48] and Messiter [49] determined the triple-deck structure about the trailing edge point. Hence, the present paper puts the coating flow problem in a similar position to that for the trailing edge problem after Goldstein's development but before that of Stewartson and Messiter: a much more complex asymptotic structure near the blade or die is needed if the analysis is to be complete or extendable to higher order. Extension to either higher-order terms (in the present expansions) or larger parameter values may need a much more complicated interaction region. The idea is then is to seek a uniformly valid solution across the singularity, and thus provide the flow details very close to the edge. In analogy to the flow near the trailing edge, we anticipate the existence of a very small region near the edge of the blade or the die where derivatives of the flow variables are of the same order in both directions. In other words, both shearing and elongation are dominant mechanisms. In this case, the correct approximate (and not just the boundary-layer) solution of the full Navier–Stokes equations must be sought. Different attempts were made early on to achieve uniform validity but the approach developed by Stewartson [48] and Messiter [49], or the now well-established triple-deck approach, looks the most promising for future development.

**Author Contributions:** This article is based on preliminary work done by M.T.H. for his MASc thesis under the supervision of R.E.K. The work was then completely rewritten and calculations carried out by R.E.K. Conceptualization, R.E.K.; methodology, R.E.K.; software, R.E.K.; validation, M.T.H.; formal analysis, R.E.K. and M.T.H.; writing—original draft preparation, M.T.H.; writing—review and editing, R.E.K.; supervision, R.E.K.; project administration, R.E.K.; funding acquisition, R.E.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** We would like to acknowledge the funding provided by the Natural Sciences and Engineering Council of Canada.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Matching between the Free-Surface and Core Layers**

The matching rule employed by Van Dyke [57] is adopted here for the streamfunction and the pressure, but the pressure will turn out to be uniformly valid across the core and the free-surface layer. Thus, only the matching of the streamfunction is detailed here. In this case, Van Dyke's rule states that

$$\mathbf{E\_nH\_m\psi} = \mathbf{H\_mE\_n\psi} \tag{A1}$$

where m and n are integers. Here, E<sup>n</sup> is the core-expansion operator, which truncates immediately after the term of order ε <sup>n</sup> where the expansion is expressed in terms of core variables. H<sup>m</sup> is the corresponding boundary-layer expansion operator. The left- and right-hand sides of expressions (A1) must be exactly the same for all the values of m and n. Depending on the values of m and n, we need different levels of matching to obtain the boundary conditions for the inner and core solutions, and simultaneously to determine the free surface height to each order in ε.

We recall from (1) that, to leading order, the stream function in the core region is ψ<sup>0</sup> (z), which can be expressed in terms of y = z − ζ(x). In this case, the stream function expression (52) for the core flow must be cast in terms of y, and takes the following form (Khayat 2014, Tillett 1968):

$$
\psi(\mathbf{x}, \mathbf{z}) = \left(\frac{1-\mathbf{G}}{2}\right) \mathbf{(y} + \varepsilon \mathbf{h}\right)^2 + \frac{\mathbf{G}}{3} (\mathbf{y} + \varepsilon \mathbf{h})^3 + \varepsilon^3 \overline{\psi}(\mathbf{x}, \mathbf{y} + \varepsilon \mathbf{h}) + \mathcal{O}\left(\varepsilon^3\right). \tag{A2}
$$

A useful expression is obtained upon expanding ψ(x, z) = ψ(x, y + εh) about y = 0:

$$\begin{array}{ll} \overline{\Psi}(\mathbf{x}, \mathbf{y} + \varepsilon \mathbf{h}) = & \overline{\Psi}(\mathbf{x}, \mathbf{0}) + \mathbf{y} \overline{\Psi}\_{\mathbf{y}}(\mathbf{x}, \mathbf{0}) + \frac{\mathbf{y}^{2}}{2} \overline{\Psi}\_{\mathbf{yy}}(\mathbf{x}, \mathbf{0}) + \dots \\ & + \varepsilon \mathbf{h} \overline{\Psi}\_{\mathbf{y}}(\mathbf{x}, \mathbf{0}) + \varepsilon \mathbf{h} \mathbf{y} \overline{\Psi}\_{\mathbf{yy}}(\mathbf{x}, \mathbf{0}) + \dots \\ & + \frac{(\varepsilon \mathbf{h})^{2}}{2} \overline{\Psi}\_{\mathbf{yy}}(\mathbf{x}, \mathbf{0}) + \dots \end{array} \tag{A3}$$

We first consider the matching between the free-surface and core layers to O(ε). The asymptotic form (27)–(29) for Ψ<sup>2</sup> is determined by considering the application of (A1) for m = 2 and n = 0. Applying E<sup>0</sup> on (A2) gives

$$\mathbf{E}\_0 \boldsymbol{\upmu} = \left(\frac{1-\mathbf{G}}{2}\right) \mathbf{y}^2 - \frac{\mathbf{G}}{3} \mathbf{y}^3. \tag{A4}$$

As this expression must be cast in inner variables when the operator H<sup>2</sup> is applied, it is rewritten in the following form:

$$\mathbf{E}\_0 \boldsymbol{\upmu} = (\mathbf{1} - \mathbf{G}) \frac{\varepsilon^2 \boldsymbol{\upmu}^2}{2} + \mathbf{G} \frac{\varepsilon^3 \boldsymbol{\upmu}^3}{3}. \tag{A5}$$

Therefore,

$$\mathrm{H}\_{2}\mathrm{E}\_{0}\psi = \left(\frac{1-\mathrm{G}}{2}\right)\mathrm{e}^{2}\eta^{2} = \left(\frac{1-\mathrm{G}}{2}\right)\mathrm{p}^{2}.\tag{A6}$$

On the other hand, to leading order, the expansion for the stream function in the free-surface layer is obtained from (20) as ψ = ε <sup>2</sup>Ψ2. Thus, E0H2ψ = ε <sup>2</sup>Ψ2, which, when matched with (A5), leads to Ψ<sup>2</sup> ∼ (1 − G) η 2 2 for large η or condition (27)–(29).

Next, we determine h0(x) in (19) by considering the application of (A1) for m = 2 and n = 1. Applying E<sup>1</sup> on (A2), we have

$$\mathbf{E}\_1 \boldsymbol{\Psi} = \left(\frac{1-\mathbf{G}}{2}\right) \mathbf{\hat{y}}^2 + 2\epsilon \mathbf{h}\_0 \mathbf{y} \mathbf{\hat{y}} + \frac{\mathbf{G}}{3} \mathbf{\hat{y}}^3 + 3\epsilon \mathbf{h}\_0 \mathbf{y}^2 \mathbf{\hat{y}}.\tag{A7}$$

The terms surviving to O ε 2 are identified by expressing E1ψ in terms of η = y ε and using (A2) for m = 1 to yield:

$$\mathbf{H}\_{2}\mathbf{E}\_{1}\boldsymbol{\psi} = \left(\frac{\mathbf{1}-\mathbf{G}}{2}\right)\mathbf{y}^{2} + \varepsilon(\mathbf{1}-\mathbf{G})\mathbf{h}\_{0}\mathbf{y}.\tag{A8}$$

On the other hand, we note from (29) that H2ψ may be written in terms of the core variables as

$$\mathrm{H}\_{2}\psi = \varepsilon^{2}\xi^{2/3}\mathrm{f}\_{2} \sim \left(\frac{1-\mathrm{G}}{2}\right)\varepsilon^{2}\xi^{2/3}(\mathfrak{G}+\mathfrak{c}\_{1})^{2} = \left(\frac{1-\mathrm{G}}{2}\right)\varepsilon^{2}\xi^{2/3}\left(\mathfrak{n}\xi^{-1/3}+\mathfrak{c}\_{1}\right)^{2},\tag{A9}$$

which yields

$$\mathbf{E}\_1 \mathbf{H}\_2 \boldsymbol{\psi} = \left(\frac{\mathbf{1} - \mathbf{G}}{2}\right) \mathbf{y}^2 + \varepsilon (\mathbf{1} - \mathbf{G}) \mathbf{c}\_1 \mathbf{x}^{1/3} \mathbf{y}. \tag{A10}$$

Equating (A8) and (A10) leads to

$$\mathbf{h\_0(x) = c\_1 x^{1/3}} = \frac{\mathbf{d\_1}}{\left(1 - \mathbf{G}\right)^{1/3}} \mathbf{x}^{1/3}.\tag{A11}$$

Even to this order, we begin to sense the difference in behavior of the film height (thickness) between slot and blade coating for G > 0 and G < 0, respectively. Clearly, (A11) indicates that in both cases the height grows like x 1/3 but at a higher rate, leading ultimately to a thinner film for slot than for blade coating. A more accurate prediction of the film thickness will be achieved by considering the next order.

In order to determine the asymptotic behavior of Ψ3(ξ, η → 0) in the free-surface layer, we set n = 0 and m = 3, and match the expressions:

$$\mathbf{H}\_3 \mathbf{E}\_0 \boldsymbol{\psi} = 2\mathbf{y}^2 - \frac{4}{3} \mathbf{y}^3,\tag{A12}$$

$$\mathbf{E}\_0 \mathbf{H}\_3 \boldsymbol{\upmu} = \mathbf{E}\_0 \left( \varepsilon^2 \mathbf{\varPsi}\_2 + \varepsilon^3 \mathbf{\varPsi}\_3 \right) \tag{A13}$$

yielding Ψ<sup>3</sup> ∼ G 3 η <sup>3</sup> or condition (35).

Finally, to obtain h1(x) and ψ(x, 0), (A1) is applied for m = n = 3. This step is algebraically much more involved, and is only summarized here. Noting that H3ψ = ε <sup>2</sup>Ψ<sup>2</sup> + ε <sup>3</sup>Ψ<sup>3</sup> = ε 2ξ 2/3 f<sup>2</sup> + ε <sup>3</sup>ξf<sup>3</sup> and using the asymptotic forms (29) and (32) give

$$\begin{split} \mathsf{H}\_{3}\mathsf{\upmu} &= \varepsilon^{2} \xi^{2/3} \Biggl\{ \big{(} 1 - \mathsf{G} ) \frac{\left( \mathfrak{n} \xi^{-1/3} + \mathsf{c}\_{1} \right)^{2}}{2} \Biggr\} \\ &+ \varepsilon^{3} \mathfrak{k} \Bigl[ \frac{\mathsf{G}}{3} \Bigl( \left( \mathfrak{n} \xi^{-1/3} + \mathsf{c}\_{1} \right)^{3} - \frac{6}{1 - \mathsf{G}} \Biggr) + \mathsf{c}\_{2} \Bigl( \mathfrak{n} \xi^{-1/3} + \mathsf{c}\_{1} \Bigr) \Bigr]. \end{split} \tag{A14}$$

Thus,

$$\begin{array}{l} \mathrm{E\_3H\_3\psi} \quad = (\mathrm{I} - \mathrm{G})\frac{\mathrm{y}^2}{2} + \frac{\mathrm{C}}{3}\mathrm{y}^3 + \varepsilon\mathrm{x}^{1/3}\mathrm{yc}\_1[(\mathrm{I} - \mathrm{G}) + \mathrm{Gy}] \\ \quad + \varepsilon^2\mathrm{x}^{2/3}\left[\frac{(\mathrm{I} - \mathrm{G})\mathrm{c}\_1^2}{2} + \mathrm{Gy}\mathrm{c}\_1^2 + \mathrm{c}\_2\mathrm{y}\right] + \varepsilon^3\mathrm{x}\left[\mathrm{c}\_2\mathrm{c}\_1 + \frac{1}{3}\mathrm{\{G}\mathrm{c}\_1^3 - \frac{\mathrm{f}\mathrm{G}}{1 - \mathrm{G}}\}\right]. \end{array} \tag{A15}$$

Recalling (A2) and (A3), we observe that

$$\begin{array}{ll} \mathrm{E}\_{3}\boldsymbol{\upmu} &= \left(\frac{\mathrm{1}-\mathrm{G}}{2}\right) \Big[ \mathrm{y}^{2} + 2\varepsilon \Big( \mathrm{h}\_{0} + \varepsilon \mathrm{h}\_{1} + \varepsilon^{2} \mathrm{h}\_{2} \Big) \mathrm{y} + \varepsilon^{2} \Big( \mathrm{h}\_{0}^{2} + 2\varepsilon \mathrm{h}\_{0} \mathrm{h}\_{1} \Big) \Big] \\ &+ \frac{\mathrm{G}}{3} \Big( \mathrm{y}^{3} + 3\varepsilon \Big( \mathrm{h}\_{0} + \varepsilon \mathrm{h}\_{1} + \varepsilon^{2} \mathrm{h}\_{2} \Big) \mathrm{y}^{2} + 3\varepsilon^{2} \Big( \mathrm{h}\_{0}^{2} + 2\varepsilon \mathrm{h}\_{0} \mathrm{h}\_{1} \Big) \mathrm{y} + \varepsilon^{3} \mathrm{h}\_{0}^{3} \Big) + \varepsilon^{3} \overline{\psi}(\mathrm{x}, 0). \end{array} \tag{A16}$$

Consequently,

$$\begin{split} \mathbf{H}\_{3}\mathbf{E}\_{3}\boldsymbol{\upmu} &= (\mathbf{1}-\mathbf{G})\frac{\mathbf{y}^{2}}{2} + \mathbf{G}\frac{\mathbf{y}^{3}}{3} + \varepsilon\mathbf{y}(\mathbf{1}-\mathbf{G}+\mathbf{G}\mathbf{y})\mathbf{h}\_{0} \\ &+ \varepsilon^{2} \Big[ \Big(\frac{\mathbf{1}-\mathbf{G}}{2} + \mathbf{G}\mathbf{y}\Big) \mathbf{h}\_{0}^{2} + (\mathbf{1}-\mathbf{G})\mathbf{y}\mathbf{h}\_{1} \Big] + \varepsilon^{3} \Big[ (\mathbf{1}-\mathbf{G})\mathbf{h}\_{0}\mathbf{h}\_{1} + \mathbf{G}\frac{\mathbf{h}\_{0}^{3}}{3} + \overline{\Psi}(\mathbf{x},\mathbf{0}) \Big]. \end{split} \tag{A17}$$

Equating (A15) and (A17), and recalling that h<sup>0</sup> = c1x 1/3 , the correction for the free surface height to the next order is obtained as

$$\mathbf{h\_1(x)} = \frac{\mathbf{c\_2}}{(1 - \mathbf{G})} \mathbf{x}^{2/3}. \tag{A18}$$

In addition, one has

$$
\overline{\psi}(\mathbf{x},0) = -\frac{\mathbf{2G}}{1-\mathbf{G}}\mathbf{x}\tag{A19}
$$

Condition (A19) yields the third boundary condition in (59)–(61).

**Appendix B. Values of Eigenvalues and Coe**ffi**cients for the First Six Modes**


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
