2.1 Dimensionless numbers

The fundamental problem of a confined dip coating system is determining the dependence of the film thickness $h_{\infty }^{\ast }$ on the system parameters ( $l$ , $u_{s}$ and $g$ ), and the physical properties of the fluid ( $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D70E}$ ). In the present study, it is assumed that a substrate is a flat plate and the coating liquid is a pure Newtonian fluid (i.e. there are no surfactants or particles) that wets the plate completely. For the flow where the foregoing assumptions are valid, the Buckingham $\unicode[STIX]{x03C0}$ theorem yields four dimensionless numbers describing the flow system uniquely. The numbers are the dimensionless thickness $T$ , the dimensionless channel width $L$ , the capillary number $Ca$ and the material number $m$ :

Note that three dimensionless numbers are required to define an unconfined dip coating flow uniquely (Tallmadge & Soroka Reference Tallmadge and Soroka1969); however, the newly introduced parameter, the channel width $l$ , generates the additional dimensionless number $L$ . Another important number, the Reynolds number $Re$ , can be obtained by the combination of the other numbers (2.2):

Here, we limit our study to a Stokes flow ( $Re\ll 1$ ); $T$ is solely determined by $Ca$ and $L$ . Accordingly, we set the value of $m$ to $10^{-5}$ in FE computations. Under this condition, we can assume a confined dip coating flow to be a Stokes flow because we explore the range in $Ca$ up to $10^{2}$ such that $Re\simeq mCa$ .

2.2 Governing equations and boundary conditions

We consider a steady flow in a confined dip coating system where the flat plate moves vertically with a constant velocity $\boldsymbol{v}_{s}^{\ast }$ (the asterisk denotes the corresponding dimensional variables throughout the present paper), as shown in figure 2, and only the part of the meniscus closer to the moving plate is examined in the present study. When the inertia is negligible, the flow is governed by the Stokes equations:

where $\boldsymbol{v}^{\ast }$ is the velocity, $\unicode[STIX]{x1D64F}^{\ast }$ is the total stress tensor and $\boldsymbol{g}$ is the gravitational field. The no-slip conditions at the moving plate and stationary wall are

The stress balance and the kinematic condition at the meniscus are

where $\boldsymbol{n}_{m}$ is the outward unit normal vector and is given by

where $\boldsymbol{i}$ , $\boldsymbol{j}$ are the unit vectors in the positive directions of the axes; and $\unicode[STIX]{x1D705}_{}^{\ast }$ is the curvature of meniscus:

Here, we set the pressure on the gas side to zero. Finally, the stress tensor $\unicode[STIX]{x1D64F}^{\ast }$ of a pure Newtonian fluid is given by

where $p^{\ast }$ is the pressure and $\unicode[STIX]{x1D644}$ is the identity tensor.

Figure 2.

Formulation of the problem. All variables are scaled for the static meniscus region. A meniscus profile $h$ is represented as a function of the $x$ coordinate. The dotted curve represents an imaginary static meniscus with a contact angle $\unicode[STIX]{x1D703}_{c}$ in the absence of flow. The solid line represents the meniscus shape due to the flow. The meniscus shape deviates from the shape of the static meniscus in the dynamic meniscus region (region II). In the inset, $\unicode[STIX]{x1D705}_{0}$ is the curvature at the centre of the imaginary static meniscus and is discussed in § 3.

2.3 Dynamic meniscus region

In this region, the balance among the viscous drag, surface tension and gravity is important, and the scaling for variables need to be chosen appropriately to examine the balance. Following Park & Homsy (Reference Park and Homsy1984), we non-dimensionalise the Stokes equations and boundary conditions by the following scales:

The variables in these scales are denoted by overbars, and the equations and conditions in the component form become

The importance of gravity in the dynamic meniscus region is represented by

which appears in the first equation of (2.15).

2.4 Static meniscus region

In this region, when the viscous force is negligible, the meniscus shape is determined by the balance between gravity and surface tension. Again, following Park & Homsy (Reference Park and Homsy1984), we use the following scales:

The variables in these scales are denoted by dropping the asterisks and the Stokes equations and boundary conditions are non-dimensionalised:

2.5 Channel region

In a confined pool, the flow inside the parallel channel between the moving plate and the stationary wall becomes significant especially during fast withdrawal. In this channel region, the pressure gradient is mainly balanced by the viscous drag. Therefore, the characteristic viscous stress is an appropriate scale for the pressure, and the remaining variables have the same scales for the static meniscus region:

The Stokes equations and boundary conditions are then non-dimensionalised as

The subscript $c$ is added to $p$ to distinguish it from that for the static meniscus region. Similar to the dynamic meniscus region, the term associated with gravity appears in the momentum balance equation (2.28):

3.1 Simple model for meniscus-controlled regime

In the low- $Ca$ regime ( $Ca^{1/3}\ll 1$ ), the entrained film is extremely thin, and hence $l_{d}\sim (h_{\infty }^{\ast }l)^{1/2}$ is small. Therefore, the meniscus remains static except for the dynamic meniscus region.

In this regime, the gravitational drainage $G_{d}$ is negligible, and the velocity field can be determined by solving (2.14)–(2.17):

Knowing the velocity field, from the kinematic condition (2.18), the mass conservation for the dynamic meniscus region can be obtained:

The term inside the parentheses of the preceding equation is the flow rate through the dynamic meniscus region. This flow rate is constant with respect to $\overline{x}$ and should be equal to that through the film entrainment region $\overline{h}_{\infty }$ . This flow-rate-matching condition yields a third-order ordinary differential equation of the meniscus profile in the dynamic meniscus region:

In addition, $\overline{h}$ in the limit to the film entrainment region is the film thickness $\overline{h}_{\infty }$ :

It becomes the boundary condition for (3.5), and the equation can be transformed into the universal form (Bretherton Reference Bretherton1961). The second derivative of the meniscus profile $\overline{h}$ in the limit to the static pool region can be obtained via numerical integration of the transformed equation:

Under the same condition ( $Ca^{1/3}\ll 1$ ), the viscous force is negligible in the static meniscus region. Therefore, the Stokes equations and boundary conditions (2.22)–(2.25) for that region can be greatly simplified as

and the integration of the foregoing equation yields

where $C$ is an integration constant. Again, we consider only the part of the meniscus along which $\text{d}^{2}h/\text{d}x^{2}>0$ . For an unconfined pool, $C$ vanishes in our coordinate system. However, for a confined pool, $C$ has a finite value because the dimensionless curvature of a static meniscus does not vanish except for $\unicode[STIX]{x1D703}_{c}=\unicode[STIX]{x03C0}/2$ , as shown in figure 2. $C$ is determined by the pressure profile along the meniscus. The liquid-side pressure $p_{0}$ at the centre point $\widetilde{O}$ of the static meniscus is selected as the pressure datum and is equivalent to $-\unicode[STIX]{x1D705}_{0}$ , where $\unicode[STIX]{x1D705}_{0}$ is the curvature at point $\widetilde{O}$ . The liquid-side pressure at any point $A$ on the meniscus in the inset of figure 2 can be obtained by considering either the hydrostatic pressure at $A$ or the capillary pressure jump across the curved meniscus at that point. Both should be identical, and hence $C=\unicode[STIX]{x1D705}_{0}$ . Equation (3.9) then becomes

Following the idea of Landau & Levich (Reference Landau and Levich1942), we also use the classical curvature-matching condition to determine the dimensionless film thickness $T$ , i.e. the curvature of the static meniscus at the ‘upper’ boundary is matched to that of the dynamic meniscus at the ‘lower’ boundary. For the static meniscus, the slope of the meniscus vanishes at the matching point:

This gives the curvature at that point:

For the dynamic meniscus, equation (3.7) can be rewritten using the variables for the static meniscus region:

Finally, from the matching between the curvatures (3.12) and (3.13), $T$ is given by

We define the term involving the curvature $\unicode[STIX]{x1D705}_{0}$ in the preceding equation as

Using $C_{c}$ , the dimensional film thickness $h_{\infty }^{\ast }$ is expressed as

The expressions inside the parentheses are identical to (2.1) with the corresponding characteristic lengths. This identity implies that the previous result of Landau & Levich (Reference Landau and Levich1942) can be corrected by introducing $C_{c}$ for confined dip coating. Therefore, we call $C_{c}$ the confinement-correction factor for the meniscus-controlled regime.

The value of $\unicode[STIX]{x1D705}_{0}$ is required to determine $C_{c}$ and the complete equation (3.14). Equation (3.10) can be integrated and rearranged as

The integration of (3.17) from $\widetilde{O}$ to $S$ (see figure 2) yields the implicit equation,

where

and $F$ and $E$ are the elliptic integrals of the first and second kind, respectively. For a given $L$ , $k$ can be evaluated from (3.18), and then we can obtain $C_{c}$ directly instead of $\unicode[STIX]{x1D705}_{0}$ .

Figure 3 shows the confinement-correction factor $C_{c}$ as a function of the dimensionless channel width $L$ for different contact angles $\unicode[STIX]{x1D703}_{c}$ , and $C_{c}$ exhibits small- and large- $L$ limiting behaviours clearly. In the small- $L$ limit, a more confined pool (smaller $L$ ) yields a meniscus with a smaller $C_{c}$ at a given $\unicode[STIX]{x1D703}_{c}$ . In the large- $L$ limit, $C_{c}$ approaches unity, which implies that confinement effects disappear, i.e. equation (2.1) is recovered from (3.16). A detailed examination of these two limits can be found in appendix C.

Figure 3.

Confinement-correction factor $C_{c}$ as a function of dimensionless channel width $L$ . The values are computed from (3.18). The arrow denotes the direction of increasing $\unicode[STIX]{x1D703}_{c}$ . The value of $C_{c}$ decreases as the liquid wets the container wall, i.e. for lower $\unicode[STIX]{x1D703}_{c}$ .

Interestingly, the power-law dependence of $C_{c}$ on $L$ is independent of $\unicode[STIX]{x1D703}_{c}$ at small $L$ , and hence $\unicode[STIX]{x1D705}_{0}$ depends only on $\unicode[STIX]{x1D703}_{c}$ . It is related to the influence of gravity on the meniscus shape and is discussed in § 6.

It is particularly interesting to analyse (3.14) in the small- $L$ limit. In this limit, $C_{c}\rightarrow L/\sqrt{2}\cos \unicode[STIX]{x1D703}_{c}$ , as explained in § C.1. The resulting equation becomes independent of $L$ and exhibits the power-law behaviour $O(Ca^{2/3})$ :

Therefore, the prefactor of (3.2) depends only on $\unicode[STIX]{x1D703}_{c}$ in the small- $L$ limit. When the coating liquid wets the wall almost completely (small $\unicode[STIX]{x1D703}_{c}$ ), the meniscus is distorted further, and a thinner film is entrained. Similar results can be found in the inclined unconfined dip coating problem (Wilson Reference Wilson1982); the thinner film is located on the side where the meniscus shape is concave.

The matching condition of Landau & Levich (Reference Landau and Levich1942) that yields (3.14) appears to be irrational. However, matched asymptotic expansions show that this condition is valid for $Ca^{1/3}\ll 1$ (Wilson Reference Wilson1982; Park & Homsy Reference Park and Homsy1984; Münch Reference Münch2002). Equation (3.14) is equivalent to the leading-order approximation of the asymptotic expansion in $Ca^{1/3}$ . Hence, the error behaviour of the equation can be predicted by investigating the previous study of Wilson (Reference Wilson1982), in which unconfined dip coating is discussed. In his study, when the inclination angle $\unicode[STIX]{x1D6FC}=0$ (vertical withdrawal), $C_{2}$ in (30) corresponds to our $C$ in (3.9) and (3.14) becomes the leading-order term of his equation (39). In terms of our characteristic length, the first correction is $O(Ca)$ (note that the first correction is $O(Ca^{1/2})$ in his study because of the different characteristic length from the present study). Therefore, one might expect the error behaviour of (3.14) to be $O(Ca)$ for low $Ca$ . However, Wilson (Reference Wilson1982) also discussed the non-uniformity in his asymptotic series that arises from $\unicode[STIX]{x1D6FC}$ . We, therefore, can predict the non-uniformity arising from the contact angle $\unicode[STIX]{x1D703}_{c}$ , which changes the meniscus profile at the given $Ca$ and $L$ , in confined dip coating. The error behaviour is briefly discussed in § 5.3.

The approach that leads to (3.14) can also be easily extended to other film-thickness-predicting equations on the basis of the curvature matching for unconfined dip coating such as that of White & Tallmadge (Reference White and Tallmadge1965) and Spiers et al. (Reference Spiers, Subbaraman and Wilkinson1974). See appendix B for details.

4.1 Limiting behaviour of film thickness in channel-controlled regime

The flow in the channel region is fully developed, as mentioned in § 2, and hence the Stokes equations with boundary conditions (2.27)–(2.30) can be greatly simplified to give a velocity field

We can obtain the flow rate in units of $u_{s}\,l$ through the channel from the velocity field:

The flow rate $q_{c}$ should be equal to that through the film entrainment region. It has been given by Van Rossum (Reference Van Rossum1958) as

The term involving $G_{c}$ represents the gravitational drainage in the film entrainment region. The drainage is negligible at high $Ca$ $\gg L^{2}$ because $G_{c}T^{2}<L^{2}/Ca$ in confined dip coating, where the film thickness is expected to be smaller than the channel width; i.e. $T<1$ . Therefore, the dimensionless thickness $T$ is given by

This result means that $T$ can be determined by considering only the pressure gradient inside the channel.

Besides, the flow inside the channel at high $Ca$ resembles the coating flow under a knife or blade. The well-developed theory for knife coatings shows that the film thickness is half the coating gap for a blade parallel to a moving plate under negligible surface tension (Ruschak Reference Ruschak1985; Coyle Reference Coyle1997). Here, the gap is the distance between the moving plate and the blade face, which is essentially identical to the channel width $l$ in confined dip coating. We can therefore expect that $\text{d}p_{c}/\text{d}x$ will eventually vanish in the high- $Ca$ limit, from knife coating theory. As a result, the flow becomes an almost pure Couette flow, and $T$ tends to $1/2$ when $Ca$ approaches infinity:

This result appears to be similar to that of unconfined dip coating: $h_{\infty }^{\ast }/l_{f}$ becomes independent of $Ca$ in the high- $Ca$ limit in the absence of inertia and surfactants (Esmail & Hummel Reference Esmail and Hummel1975b ; Weinstein & Ruschak Reference Weinstein and Ruschak2001; Jin et al. Reference Jin, Acrivos and Münch2005). However, the film thickness $h_{\infty }^{\ast }$ itself increases without bound with $Ca$ in unconfined dip coating, whereas $h_{\infty }^{\ast }$ is bounded in confined dip coating.

5.1 Finite element model

To solve the Stokes equations with the Galerkin finite element method (G/FEM), we need more boundary conditions, in addition to conditions (2.6)–(2.9). These boundary conditions are represented in the following.

A stationary wall with the no-slip condition should be included in the analysis to incorporate the confinement effects. Such a wall can be found in early numerical models for unconfined dip coating (Jenny & Souhar Reference Jenny and Souhar2009; Filali et al. Reference Filali, Khezzar and Mitsoulis2013); however, there is also an inlet feed port to compensate for the mass loss at the outflow boundary in those models. The reliability of a numerical model with such an inlet feed would depend on the location of the feed; the flow through a feed that is not far enough from the meniscus affects the film thickness when a pool is confined. Besides, the numerical models with sufficient distance between the meniscus and the feed may be more computationally expensive than our FE model explained in the following.

The meniscus contacts the stationary wall at an angle $\unicode[STIX]{x1D703}_{}$ on the contact line:

where $\boldsymbol{n}_{m}$ and $\boldsymbol{n}_{w}$ are the outward normal unit vectors on the meniscus and the wall, respectively. The definitions of the normal and tangent vectors and a contact angle in our FE model are shown in figure 4. When the coating liquid near the stationary wall remains in static equilibrium, the static contact angle $\unicode[STIX]{x1D703}_{c}$ is used instead of $\unicode[STIX]{x1D703}_{}$ .

Figure 4.

Numerical flow domain and corresponding boundary conditions. All boundary conditions, unit vectors and a contact angle are represented. The free boundary condition is applied to both the inflow (bottom) and outflow (top) boundaries. The velocity profile is not prescribed.

Dip coating is a self-metered method in that the flow rate is determined by the force balance inside the coating flow (Weinstein & Ruschak Reference Weinstein and Ruschak2001). Therefore, a velocity profile should not be prescribed at the inflow or outflow boundary a priori; such a prescribed velocity profile presets the film thickness $h_{\infty }^{\ast }$ .

In this regard, the FEM has an advantage over other methods, because of the free boundary condition (Papanastasiou, Malamataris & Ellwood Reference Papanastasiou, Malamataris and Ellwood1992). The FEM can extend the validity of the weak form of the governing equations to synthetic boundaries (inflow and outflow boundaries), instead of replacing them with unknown essential or natural boundary conditions, such that

where $\boldsymbol{n}_{b}$ is the outward unit normal vector at either the inflow or outflow boundary. The free boundary condition (5.2), which is typically used for an outflow boundary, appears not to impose any boundary condition at all; however, it provides an effective boundary condition (Renardy Reference Renardy1997).

We apply the free boundary condition (5.2) not only to the outflow boundary but also to the inflow boundary. To the best of the authors’ knowledge, this has never been attempted before. In our FE model, both the inflow and outflow boundaries are located sufficiently far from each other so that the velocity fields and the gradient of the pressure fields are expected to approach asymptotic limits at both ends. This asymptotic state includes ‘fully developed’ flows (Renardy Reference Renardy1997). Applying the free boundary condition to the inflow boundary allows the pressure to take any value at the boundary. It yields the undetermined position of the contact line with the contact angle prescribed by (5.1). However, at a given $Ca$ and $L$ , both the position of the contact line and the contact angle should be fixed to have a unique solution in the frame of reference for our model, but it is impossible to impose them as boundary conditions in our FE framework simultaneously. Therefore, an additional condition or constraint is required to find the unique solution. It would seem that the mesh configuration of the physical domain of interest provides the additional condition required to determine the position of the contact line implicitly. When the mesh generation equation (explained in the following) is solved with (2.4) and (2.5) simultaneously, we can find the unique (steady-state) solution. However, the detailed mechanism of how the mesh configuration fixes the position of the contact line is not clear.

As discussed earlier, the physical domain of interest of the confined dip coating flow should be surrounded by three physical boundaries (the meniscus, the stationary wall and the moving plate) and two synthetic boundaries (the inflow and outflow planes). All the boundary conditions are summarised in figure 4.

To find the meniscus shape via standard techniques for boundary value problems, the set of differential equations and boundary conditions on the unknown, deformable ‘physical’ region $R$ in the $x^{\ast }y^{\ast }$ -plane should be transformed one to one into an equivalent set on the known, fixed ‘computational’ region $G$ in the $\unicode[STIX]{x1D709}\unicode[STIX]{x1D702}$ -plane. Here $R$ and the set on it are transformed by an equation of the form $\boldsymbol{x}^{\ast }=\boldsymbol{M}(\unicode[STIX]{x1D743})$ , where $\boldsymbol{x}^{\ast }=(x^{\ast },y^{\ast })$ and $\unicode[STIX]{x1D743}=(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ . Any function $f(\boldsymbol{x}^{\ast })$ defined on $R$ can be thought of as a function $f(\boldsymbol{x}^{\ast }(\unicode[STIX]{x1D743}))$ defined on $G$ as well. The image $\unicode[STIX]{x1D743}$ of $\boldsymbol{x}^{\ast }$ under the inverse transformation $\boldsymbol{M}^{-1}$ is a solution to elliptic partial differential equations, which is similar to that for the diffusional mass transport for the dilute regime,

where $D_{\unicode[STIX]{x1D709}}$ and $D_{\unicode[STIX]{x1D702}}$ are mesh diffusivities, which determine the gradients of the coordinate potentials $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D702}$ adjusting the nodal distribution inside the ‘physical’ region $R$ . The nodal distribution at the boundaries, prescribed by the functions defining their geometries, is specified by a stretching function (Vinokur Reference Vinokur1983). The iso-contour curves of $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D702}$ define the boundaries of the quadrilateral elements used to tessellate $R$ , and the intersection points of these curves determine the nodal positions. The meniscus shape is implicitly determined by the kinematic condition (2.8). Equations (5.3) have been called the elliptic mesh generation equations, and they have been discussed in detail by de Santos (Reference de Santos1991).

5.2 Solution using the finite element model

The field variables, such as the velocity, pressure and nodal position, are approximated by a linear combination of the polynomial basis functions. The velocity and nodal position are approximated by a Lagrangian biquadratic function $\unicode[STIX]{x1D719}_{i}^{}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ , and the pressure is approximated by a linear discontinuous function $\unicode[STIX]{x1D713}_{k}^{}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ :

where $n$ is the number of nodes used for the nodal position and velocity and $m$ is the number of nodes used for the pressure, and the coefficients $\boldsymbol{x}_{i}^{\ast }$ , $\boldsymbol{v}_{i}^{\ast }$ and $p_{i}^{\ast }$ are unknown values to be determined at each node.

The weak forms of (2.4), (2.5) and (5.3) are obtained by multiplying them by the weighting functions and then integrating over the ‘physical’ region $R$ . The essential boundary conditions are implemented by replacing the corresponding weighted residual equation with the desired velocity or node specification. The natural boundary conditions are implemented through the boundary integrals, which follow from the divergence theorem. Note that the free boundary condition (5.2) is treated as the natural boundary condition.

The G/FEM reduces the Stokes and mesh generation equations (2.4), (2.5) and (5.3) to a set of nonlinear algebraic equations with respect to the basis function coefficients,

where $\boldsymbol{z}$ is the solution vector consisting of the FE coefficients of $\boldsymbol{x}_{i}^{\ast }$ , $\boldsymbol{v}_{i}^{\ast }$ and $p_{i}^{\ast }$ ; and $\unicode[STIX]{x1D740}$ is the parameter vector containing the parameters for the dip coating flow, such as $l$ , $u_{s}$ , $\unicode[STIX]{x1D703}_{c}$ , $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D70E}$ . The system of nonlinear equations (5.5) is solved using Newton’s method:

where $\unicode[STIX]{x1D645}_{\boldsymbol{R}}(\boldsymbol{z})$ is the Jacobian matrix of $\boldsymbol{R}$ . The solution vector $\boldsymbol{z}$ is accepted as the solution to the problem when $\Vert \boldsymbol{R}(\boldsymbol{z}^{i})\Vert _{2}<5\times 10^{-9}$ . The typical value of the solution update $\Vert \unicode[STIX]{x1D6FF}\boldsymbol{z}^{i}\Vert _{2}$ is of the order of $10^{-5}$ at convergence.

We use a direct numerical solver based on the frontal algorithm (Duff, Erisman & Reid Reference Duff, Erisman and Reid1989). It appears that there is no zero pivot during the matrix decomposition step; thus, we can safely conclude that the free boundary condition (5.2) at both the inflow and outflow boundaries does not yield an ill-posed problem in the resulting discrete-level matrix system (5.6).

The initial guess was obtained by replacing the unknown meniscus shape with a known arc of a circle and substituting the stress balance condition (2.8) with the kinematic condition (2.9) and Navier’s slip condition with a large slip coefficient $\unicode[STIX]{x1D6FD}=10^{8}$ along the meniscus, $\unicode[STIX]{x1D6FD}(\boldsymbol{t}_{m}\boldsymbol{n}_{m}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{\ast })=\boldsymbol{t}_{m}\boldsymbol{\cdot }\boldsymbol{v}^{\ast }$ . Nodal positions were not included in the unknowns because the normal and tangential unit vectors at the meniscus were known at this stage; Navier’s slip condition and the kinematic condition (2.9) were applied to the $x^{\ast }$ - and $y^{\ast }$ -components of the momentum equation (2.5), respectively. The base solution at $Ca=0.01$ , $L=5$ and $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ was computed using the initial guess. Figure 5 illustrates the mesh configuration of the base solution. Along the meniscus, a total of 125 elements and 251 nodes were used to resolve the meniscus shape.

Figure 5.

Mesh configuration at $Ca=0.01$ , $L=5$ and $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ . 2775 quadrilateral elements are tessellated into the physical domain. A total of 125 elements and 251 nodes is used to resolve the curved meniscus shape. The aspect ratio of the mesh configuration is adjusted for visualisation ( $d=6\,l$ ).

The solutions for different values of the capillary number $Ca$ and the dimensionless channel width $L$ were computed by using the pseudo-arc-length continuation (Bolstad & Keller Reference Bolstad and Keller1986), starting from the base solution. $u_{s}$ and $l$ were chosen as the continuation parameters to adjust $Ca$ and $L$ , respectively. During the continuation in $l$ , the distance $d$ between the inflow and outflow boundaries should be adjusted simultaneously to ensure a fully developed flow at the inflow and outflow boundaries. When the inflow is not fully developed, the accuracy of the flow rate through the flow domain could be compromised, and hence the film thickness $h_{\infty }^{\ast }$ . Unlike other numerical models (Jin et al. Reference Jin, Acrivos and Münch2005; Campana et al. Reference Campana, Ubal, Giavedoni and Saita2013), we did not restrict the shape of the free surface at the corner of the element where the outflow boundary and the free surface are placed. Hence, the free surface slope with respect to the moving plate is a major indicator of fully developed flow. When the outflow is not fully developed, the thickness being measured from the moving plate to the free surface along the outflow boundary is larger than the film thickness $h_{\infty }^{\ast }$ . According to our numerical tests, it turns out that $d/l=6$ ensures a virtually parabolic velocity profile and a zero slope of the free surface at the corner at the outflow boundary. We therefore maintained $d/l=6$ during the continuation.

Steady-state solutions are obtained within 3–5 iterations except for the base solution obtained using the initial guess. To guarantee a high accuracy for $h_{\infty }^{\ast }$ , we verified that the mesh configuration near the highly curved meniscus is properly resolved for every solution.

5.3 Numerical results

We validated our FE model by comparison with results from other models and experiments on unconfined dip coating. See appendix D for details.

The results of the FE computations confirm that (3.20) and (4.5) represent the meniscus- and channel-controlled regimes, respectively. The curves computed at the dimensionless channel width $L<0.5$ overlap and are virtually indistinguishable, as shown in figure 6. The overlapped curve exhibits power-law dependence and is independent of $L$ in the low- $Ca$ regime. This dependence is already predicted by (3.20).

Figure 6.

Dimensionless film thickness $T$ versus capillary number $Ca$ for $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ and various $L$ . The curves of $T$ overlap at $L<0.5$ .

Figure 7 shows the difference between the analytical and numerical results. Interestingly, the error behaviours are distinct at different contact angles. The error follows a power-law behaviour $\unicode[STIX]{x1D716}=O(Ca^{m})$ , where $m$ is unity for $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ , but changes to $m\simeq 0.75$ for $\unicode[STIX]{x1D703}_{c}=60^{\circ }$ . This distinction appears to be due to the non-uniformity mentioned in § 3.1. A full asymptotic analysis is required to investigate the non-uniformity; however, it is beyond the scope of the present study.

Figure 7.

Error $\unicode[STIX]{x1D716}$ behaviour. The error is defined as the difference between the analytical and numerical results. The dimensionless film thickness $T$ is predicted by (3.14) and (3.20) for $L=1$ and $L=0.1$ , respectively. The plot shows the errors of $O(Ca)$ at $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ and $O(Ca^{4/3})$ at $\unicode[STIX]{x1D703}_{c}=0^{\circ }$ , as indicated by Wilson (Reference Wilson1982) and Park & Homsy (Reference Park and Homsy1984), respectively.

Figure 8 shows the effect of the contact angle $\unicode[STIX]{x1D703}_{c}$ on the dimensionless thickness $T$ . $T$ increases with $\unicode[STIX]{x1D703}_{c}$ as can be seen in (3.20). However, the dependence of $T$ on $\unicode[STIX]{x1D703}_{c}$ decreases as $Ca$ increases, as shown in figure 8(a). This is because the meniscus shape becomes less sensitive to $\unicode[STIX]{x1D703}_{c}$ for high $Ca$ , as can be seen from a comparison of figures 8(b) ( $Ca=0.01$ ) and (c) ( $Ca=5$ ). The results clearly show that the capillary effect is crucial in the meniscus-controlled regime discussed in § 3, and the pressure gradient inside the channel region is crucial in the channel-controlled regime discussed in § 4.

Figure 8.

Influence of container wall wetting on film thickness and meniscus shape for $L=0.1$ and various $Ca$ . (a) Plot of $R_{T}\equiv (T-T_{0})/T_{0}$ , where $T_{0}$ is the value at $\unicode[STIX]{x1D703}_{c}=0$ . The effect of the wall wetting on the film thickness diminishes as $Ca$ increases. (b,c) The coordinates are scaled with $l$ . The arrows denote the direction of increasing $\unicode[STIX]{x1D703}_{c}$ . The meniscus shape changes significantly for $Ca=0.01$ , whereas it changes slightly for $Ca=5$ .

In the high- $Ca$ regime, the pressure gradient $\text{d}p_{c}/\text{d}x$ at the inflow boundary and the gravity $G_{c}$ monotonically decrease as $Ca$ increases, but $G_{c}$ decreases faster than $\text{d}p_{c}/\text{d}x$ ( $G_{c}=o(\text{d}p_{c}/\text{d}x)$ ), as shown in figure 9(a). Note that $\text{d}p_{c}/\text{d}x$ does not rapidly vanish as $Ca$ increases, as shown in figure 9(b). However, we can expect that $\text{d}p_{c}/\text{d}x$ will vanish eventually in the high- $Ca$ limit, and the flow will become a Couette flow (see appendix A for details). The numerical results are consistent with the discussion in § 4.1. Interestingly, $\text{d}p_{c}/\text{d}x$ is independent of $L$ at small $L$ ( ${<}0.5$ ) for high $Ca$ , as shown in figure 9(b), and hence $T$ is also independent of $L$ . This independence is identical to that found in the low- $Ca$ regime with small $L$ . Finally, all the numerical results show that $l$ is the critical parameter for both regimes.

Figure 9.

(a) The ratio of gravitational force $G_{c}$ to pressure force $\text{d}p_{c}/\text{d}x$ and (b) $\text{d}p_{c}/\text{d}x$ as a function of capillary number $Ca$ for different $L$ at $\unicode[STIX]{x1D703}_{c}=30^{\circ }$ . (a) $G_{c}$ decays faster than $\text{d}p_{c}/\text{d}x$ as $Ca$ increases. (b) When the pool is sufficiently confined ( $L<0.5$ ),  $\text{d}p_{c}/\text{d}x$ becomes independent of $L$ . Although $\text{d}p_{c}/\text{d}x$ does not vanish completely, even at an extremely high $Ca$ of up to $10^{2}$ , it clearly decreases monotonically.