Dynamics of driven three-dimensional thin films : From hydrophilic to superhydrophobic regimes

We study the forced displacement of a thin film of fluid in contact with vertical and inclined substrates of different wetting properties, that range from hydrophilic to hydrophobic, using the lattice-Boltzmann method. We study the stability and pattern formation of the contact line in the hydrophilic and superhydrophobic regimes, which correspond to wedge-shaped and nose-shaped fronts, respectively. We find that contact lines are considerably more stable for hydrophilic substrates and small inclination angles. The qualitative behavior of the front in the linear regime remains independent of the wetting properties of the substrate as a single dispersion relation describes the stability of both wedges and noses. Nonlinear patterns show a clear dependence on wetting properties and substrate inclination angle. The effect is quantified in terms of the pattern growth rate, which vanishes for the sawtooth pattern and is finite for the finger pattern. Sawtooth shaped patterns are observed for hydrophilic substrates and low inclination angles, while finger-shaped patterns arise for hydrophobic substrates and large inclination angles. Finger dynamics show a transient in which neighboring fingers interact, followed by a steady state where each finger grows independently. © 2008 American Institute of Physics. DOI: 10.1063/1.2940726


I. INTRODUCTION
The stability of forced contact lines is the key to a variety of dynamic wetting processes 1,2 that range from linear or spin coating of surfaces 3 to monodisperse drop formation 4 and in-drop mixing in microchannels. 5,6n this paper we study the dynamics of a threedimensional thin film spreading on a dry inclined substrate, where the contact line is gravity driven, as shown in Fig. 1.In such a case, the contact line destabilizes and gives rise to the growth of differently shaped patterns.
The stability of the front has been subject of active work during the last two decades.Experimentally, efforts have focused on characterizing the destabilization of the contact line and the patterns that emerge at long times.Patterns have been shown to be influenced by fluid wetting properties.Silvi and Dussan 7 already observed such a dependence, as patterns in their experiments ranged from sawtooth shaped for mostly wetting fluids to finger shaped for mostly nonwetting fluids.This feature was verified in the experimental work of Jerret and de Bruyn, 8 who found that the width of the pattern is a function of the wetting properties as well as of the imposed forcing, i.e., the inclination angle of the substrate.In general, experiments indicate that closely spaced patterns are expected for mostly wetting fluids at small inclination angles.
Theoretically, the main framework is that of the lubrication approximation of the Navier-Stokes equations, subject to some boundary condition at the moving contact line that regularizes the so-called viscous dissipation singularity. 2In the lubrication approximation one reduces the dynamics of the thin film to a two-dimensional evolution equation for the film thickness h, which is governed by a single parameter D = ͑3Ca͒ 1/3 cot͑␣͒, where Ca is the capillary number and ␣ is the inclination angle of the substrate.Wetting properties are fixed by the choice of the boundary condition at the contact line, which generally corresponds to perfectly wetting fluids, where the dynamic contact angle D is zero ͑precursor film model͒, or to mostly wetting fluids, where D is small ͑slip model͒.
The lubrication equations were first examined qualitatively by Brenner,9 who proposed that the capillary ridge located near the contact line ͑a consequence of the competition between capillary and driving forces͒ induces transverse flows that destabilize the contact line when it is perturbed.Troian et al. 10 solved numerically the linear stability problem in the long-wavelength limit using a precursor film model for the contact line.Indeed, they found that the perturbation growth rate can only be positive if the height profile is nonmonotonic.Bertozzi and Brenner 11 extended the linear stability analysis to arbitrary wavelengths using the same contact line model finding similar results.Spaid and Homsy 12 compared the stability of the front for the precursor film model and the slip model.They found that the front is rather insensitive to the contact line model used and that the front stability is very similar whenever the precursor film thickness and the slip length are comparable.Anyhow, all studies agree on the destabilization mechanism, which is based on transverse variations of the height profile in the presence of a capillary ridge.In this situation, the troughs of a perturbation feed the peaks ͑points B and A in Fig. 1, respectively͒.As a consequence, the peaks increase their local height and, in turn, their velocity.In other words, a given value of the control parameter fixes the size of the capillary ridge and, in turn, the stability of the front.
The long time evolution of the lubrication equations has been addressed numerically by Kondic and Diez 13 for completely wetting fluids and by Moyle et al. 14 for partially wetting fluids.Kondic and Diez studied the effect of the inclination angle ␣ and observed the formation of fingers at large ␣, while for sufficiently small inclination angle the growth of the pattern saturates to a sawtooth shape.Moyle et al. 14 focused on the effect of the capillary number and dynamic contact angle on the shape of fingers.For small Ca and D they observe sawtooth shaped patterns, while for larger values of both parameters they observe the growth of fingers.
In spite of these results, there are a number of situations that fall beyond the small-Ca and small-D limit of lubrication theory.
Concerning the role of large D , an interesting effect was observed by Veretennikov et al., 15 who considered the flow of films of partially wetting fluids at small inclination angles.Instead of the wedge shaped front typical of lubrication theory, they observe a multivalued film thickness profile; a nose shaped front, where the flow is no longer unidirectional, but has a vortex structure.This three-dimensional limit is typical of the hydrophobic ͑ D Ͼ 90°͒ and superhydrophobic ͑ D Ͼ 120°͒ regimes.
Alternatively to lubrication theory, three-dimensional flows of immiscible fluids can be modeled by a mesoscopic version of the Navier-Stokes equations coupled to a convection-diffusion equation for an order parameter.Such an approach has been used by a number of authors in fields as varied as spinodal decomposition 16 and shearing 17 of binary fluids, contact line motion between sheared solid plates, 18 and viscous fingering in Hele-Shaw cells. 19,20Under this approach, the contact line and the free film surface have a finite size, a feature that eliminates the complications of a sharp interface.Additionally, the evolution equation of the order parameter can be constructed from a well-known equilibrium state, thus allowing for a correct description of wetting properties.In our model, wetting is introduced by considering a surface free energy of the Cahn type, which fixes a boundary condition for the order parameter field at the solid.
In this paper we perform numerical simulations of the mesoscopic Navier-Stokes and convection-diffusion equations to study the full three-dimensional problem of a thin film flowing down a dry inclined substrate using the lattice-Boltzmann method.We focus mainly on situations in which neither Ca nor D are small, so the lubrication equations are no longer valid.We consider fluids of arbitrary wetting properties at various inclination angles.We study in detail the unperturbed front by characterizing the wedge and nose configurations.We then study the stability of these solutions where we pay particular attention to the effect of wetting properties.For long periods, we study the formation of nonlinear patterns.We find that, depending on the wetting properties and the inclination angle, saturated sawtooth patterns or steadily growing fingers are obtained, and subsequently study the robustness of the finger pattern.
The present paper is organized as follows.In Sec.II we present the three-dimensional model which consists of the coupled Navier-Stokes and convection-diffusion equations subject to a Cahn wetting boundary condition.We discuss the lubrication limit of these equations.In Secs.III and IV we present the lattice-Boltzmann method and the simulation strategy, respectively.The flat contact line problem is addressed in Sec.V, where we explore the effects of wetting, capillarity, and gravity on a propagating front.In Sec.VI we perform a numerical linear stability analysis for a slightly perturbed contact line.We explore the effect of wetting and gravity and compare to previous results. 12In Sec.VII we study the nonlinear regime of contact line dynamics.In Sec.VII A we study the growth of sawtooth patterns and fingers, which is controlled either by varying substrate wettability or substrate inclination angle.Finger growth is studied in detail in Sec.VII B, where we consider the effect of initial conditions, system size, and nonlinear interaction between neighboring fingers.In Sec.VIII we discuss our results and present the conclusions of this work.

II. GOVERNING EQUATIONS
We consider the motion of a thin film of fluid that flows down an inclined solid substrate under the action of gravity, as depicted in Fig. 1.We consider a binary fluid composed of two immiscible phases of very different viscosities.In this approximation, the more viscous fluid plays the role of the film, while the less viscous fluid plays the role of the surrounding air.Each phase is determined by a dimensionless composition variable, or order parameter, ͑r ជ͒.In equilibrium, minimizes a free energy functional given by where ͑r ជ͒ is the fluid density.The first term in the integrand, V͑ , ͒ = A 2 / 2+B 4 / 4+ ln , is a volume contribution and determines the equilibrium composition values for each phase, while the second term penalizes composition gradients by a factor , thus defining the free energy of the interface.
Minimization of F leads to the chemical potential, and total pressure tensor, where ␦ ញ is the diagonal matrix.The pressure tensor has an ideal contribution given by P ញ = ͑ / 3͒␦ ញ , and an order param- eter contribution.The equilibrium order parameter profile is * ͑n͒ =− eq tanh͑n / ͒, where eq = ͑−A / B͒ 1/2 is the bulk equilibrium value of the order parameter, = ͑− / 2A͒ 1/2 is the length scale of the interfacial region, and n is a coordinate that measures the distance from the interface in the normal direction.The energy per unit area of the interface is The divergence of the pressure tensor yields the force per unit volume that acts on the fluid: −ٌ ជ P − ٌ ជ .The first term is the pressure gradient, while the second arises from order parameter inhomogeneities.Including this last term, the Navier-Stokes equations are written as follows:

͑1͒
where v ជ is the fluid velocity, is the fluid viscosity, and g ជ is the acceleration of gravity.
The dynamics of the order parameter are described by a convection-diffusion equation, where M is a mobility.For small deviations from the equilibrium configuration, an expansion of the chemical potential in powers of − * yields a first-order diffusion coefficient D = M͑A +3B eq 2 ͒.In the sharp interface limit ͑⑀ → 0͒ Eq. ͑1͒ reduces to the usual Navier-Stokes equations in the bulk of each fluid, while at the interface the Gibbs-Thomson condition is recovered.In addition, for low D , the diffusive term in Eq. ͑2͒ vanishes.Consequently, the order parameter is advected by the velocity field.
In this limit, the standard approach of previous works 9-14 is to restrict the equations to the lubrication regime, which follows by considering only the components of the equations of motion parallel to the solid substrate averaged over the film thickness, h͑x , y͒.Consequently, the local fluid velocity is given by where the x-component of the gravitational force comes from the choice of the reference frame ͑see Fig. 1͒.In Eq. ͑3͒, the local pressure is composed of a surface term and a hydrostatic term, P =−ٌ 2 h + gh cos͑␣͒.The local velocity is then a function of the local height and of its spatial derivatives.
Averaging the continuity equation over h and combining with Eq. ͑3͒ give a fourth-order partial differential equation for the film thickness, which is the main equation in lubrication theory,

͑4͒
This equation can be nondimensionalized by using units hЈ = h / h c , ͑xЈ , yЈ͒ = ͑x / x c , y / x c ͒, and tЈ = t / t c , where h c is the film height far from the contact line.Space and time units are chosen as x c = h c ͑3Ca͒ −1/3 and t c = x c / U, where Ca= U / is the capillary number and U = h c 2 g sin͑␣͒ / 3 is the average velocity of the film.In terms of the nondimensional variables Eq. ͑4͒ becomes ‫ץ‬hЈ ‫ץ‬tЈ

ˆ͔, ͑5͒
where D = ͑3Ca͒ 1/3 cot͑␣͒ measures the magnitude of the normal gravitational force.According to lubrication theory, D is the only relevant control parameter of the system.

A. Boundary conditions
Wetting properties in the mesoscopic model presented above are included by means of an interaction energy between the fluid and the solid where S is the solid-fluid boundary and f S ͑ S ͒ is the free energy per unit surface, which is a function of the local composition S .In equilibrium this expression yields the following boundary condition at the solid-fluid interface: where n ˆis normal to the boundary.Here we consider a linear form for the surface free energy density, f S = H S , which already allows a wetting transition. 21or a rectangular box of dimensions L ϫ W ϫ b, stick boundary conditions are imposed at the solid, v ជ͑x , y , z =0͒ =0 ជ , while no flow boundary conditions are imposed for the order parameter, v ជ͑x , y , z =0͒ =0 ជ .Far above from the film and at both ends of the box the flow is homogeneous.We ensure these conditions by fixing Contact line dynamics arise from the diffuse nature of the interface, which allows for slip in the interfacial region by a diffusive mechanism.The size over which slip takes place, l D , is a function of the fluid properties and has been estimated by Briant and Yeomans, 18 who have given a scaling relation,

III. LATTICE-BOLTZMANN METHOD
We solve numerically Eqs.͑1͒ and ͑2͒ by means of the lattice-Boltzmann algorithm presented in Ref. 16.The dynamics are introduced by discretized Boltzmann equations of two velocity distribution functions, and In these equations, f i and g i are velocity distribution functions of particles that move with velocity c ជ i , where the index i counts over the model velocity set.Space is discretized as a cubic lattice where nodes are joined by velocity vectors.Space and time units in Eqs.͑7͒ and ͑8͒ are set to unity.Likewise, the density of the fluids is set to one.We use the D3Q15 velocity set, 21 which consists of 15 velocity vectors: six of magnitude 1 that correspond to nearest neighbors, eight of magnitude ͱ 3 that correspond to third-nearest neigh- bors, and one of zero magnitude that accounts for rest particles.In the D3Q15 model the speed of sound is c s =1/ ͱ 3.
In Eqs.͑7͒ and ͑8͒, distribution functions are first relaxed to equilibrium values, represented by f i eq and g i eq , with relaxation time scales f and g .The term F i f is related to the external forcing.Following the collision stage, distribution functions are propagated to neighboring sites.
Hydrodynamic variables are defined through moments of the f i and g i ; the local density and order parameter are introduced as ͚ i f i = and ͚ i g i = .The fluid momentum is defined as Local conservation of mass and momentum is enforced through the conditions ͚ i f i eq = , ͚ i g i eq = , and In equilibrium, the order parameter current, pressure tensor and chemical potential are defined as The equilibrium distribution functions and the forcing term are written as expansions in powers of v ជ, 22 i.e., Here, stands for the three possible magnitudes of the c ជ i set.
Equations ͑1͒ and ͑2͒ are recovered by performing a Chapman-Enskog expansion of Eqs.͑7͒ and ͑8͒. 23The lattice-Boltzmann scheme maps to the hydrodynamic model through the relaxation time scales, i.e., = ͑2 f −1͒ / 6 and M = ͑ g −1/ 2͒M ˆ, and through the body force f ជ = g ជ.
Solid boundaries in the lattice-Boltzmann method are implemented by means of the well-known bounce-back rules. 21,23In the fluid nodes that link to solid nodes, the propagation scheme is modified so the distribution functions are bounced back to the fluid rather than absorbed by the solid.As a consequence, a stick condition for the velocity is recovered approximately halfway from the fluid node to the solid node.

IV. SIMULATION SETUP AND MEASUREMENTS
Our aim is to study the flow of thin films for different Ca, ␣, and E .To achieve this we must choose values for the model parameters that assure the numerical stability of the algorithm.
The typical experimental situation is that of the forced spreading of a viscous fluid on a substrate in the presence of air at relatively small velocities.Inertial effects, quantified by the Reynolds number Re= h c U / , are expected to be small.This is enforced in our simulations by neglecting the nonlinear term in Eq. ͑1͒.The local viscosity is set according to the mixing rule ͑͒ = ͑ + A ͒͑1+ / eq ͒ / 2, where c = ͑ − A ͒ / ͑ + A ͒ is the viscosity contrast between the fluids and the subscript A stands for the less viscous phase.To reproduce the high asymmetry of viscous dissipation between both fluids, we fix c = 0.9, so viscous dissipation occurs mainly inside the film.
Given that we work with a finite-size interface, the ratio between the interfacial and macroscopic scales should be reasonably small so that bulk effects do not become shadowed by the interface structure.Still, the interface has to be sufficiently thick to avoid lattice artifacts.This is ensured by fixing ⑀ = 0.57 which is sufficiently large to resolve the interface. 21For this interface size, film dynamics is correctly resolved by taking h ϳ O͑10͒.
As we have mentioned before, the contact line slips by virtue of diffusion.This means that the diffusion coefficient D must be chosen to ensure that the order parameter relaxes rapidly compared to the advection time scale.In a previous study we have shown that this is accomplished as long as the product between the Péclet number Pe= Uh / D and the capillary number is O͑10 −1 ͒. 19 After these considerations all other parameters are fixed by a choice of Ca, E , and ␣.The corresponding velocity, viscosity, and surface tension lie in the ranges of 0.001ഛ U ഛ 0.01, 0.05ഛ ഛ 0.5, and 0.001ഛ ഛ 0.01, respectively.
The body force terms are fixed according to g x =3U͑ +1͒ / h 2 and g z =−g x cot͑␣͒.Then, choosing ␣ fixes the capillary number.As for wetting properties, fixing E determines the value of the wetting coefficient H.
Initial conditions are set by fixing the height h, length l, and width w of a rectangular film.Lattice nodes lying inside of the film are set to = eq , while the rest of the nodes are set to =− eq .The initial velocity of all nodes is set to ͑U ,0,0͒.
For a given configuration of the order parameter, a choice of an iso-curve determines the location of the interface.Here we choose the value = 0.The dynamic contact angle is extracted by a circle fit to interface coordinates in the xz-plane.To improve the accuracy, the fitting region is varied systematically by increasing the arc length along the interface until the fitting error becomes minimum.As is standard, we only consider the part of the interface that couples to bulk dynamics and discard the region close to the wall, where the profile distorts due to the wetting boundary condition 24 ͓see, for instance, Fig. 3͑e͔͒.

V. FLAT CONTACT LINE
Before exploring the stability and pattern formation of the thin film let us consider the effect of wetting properties on the evolution of an unperturbed front, for which the contact line is flat.Here we depart from the lubrication regime, in which the only parameter D gathers the effects of Ca and ␣.Apart from these parameters, we study fronts at three different equilibrium contact angles: E = 0°, E = 25°, and E = 90°.

A. Vertical substrates
Let us first examine the low Ca and low D limit for a vertical substrate ͑D =0͒, which corresponds to the lubrication regime, and compare to lattice-Boltzmann results.We perform a simulation fixing E = 0°and Ca= 0.04.The resulting profile, shown in Fig. 2͑a͒, has a dynamic contact angle of D = 28°.In our model, the contact line slips over the solid surface due the diffusive mechanism that originates from the finite-sized interface.Therefore, we choose to compare to the slip-velocity lubrication model of Ref. 12, which follows from coupling Eq. ͑5͒ to a slip boundary condition at the solid-fluid interface.For a flat, steady contact line, the height profile of such a model satisfies the equation where ␣ h is the slip parameter.Typical values of this parameter used in Ref. 12 lie in the range of 0.01ഛ ␣ h ഛ 0.1.Numerical solution of this equation is performed by fixing ␣ h , the contact slope C = ͑3Ca͒ −1/3 tan D , and using a fourth order Runge-Kutta integration scheme.To compare to lattice-Boltzmann results, we vary ␣ h and C until a best fit of the lattice-Boltzmann profile is achieved.Figure 2͑a͒ shows the comparison for ␣ h = 0.04 and C = 1.1.Both profiles compare acceptably.Nonetheless, in the region behind the capillary ridge a characteristic "dip" is absent from the lattice-Boltzmann solution.To improve this comparison, we decrease the diffusion strength at the contact line, which is responsible for the slip velocity and, in turn, the strength of this dip.Decreasing this strength has the effect of increasing the contact slope slightly.Nonetheless, the resulting profile now compares better to the lubrication prediction with ␣ h = 0.05 and C = 1.56, as can be seen in Fig. 2͑b͒.In both cases, deviations are comparable to the size of the interface, as shown in the figure.
Let us now explore the effect of larger equilibrium contact angles at different Ca values.Plots of the steady free surface are shown in Fig. 3. Simulation parameters for these runs are summarized in Table I.For all runs the free surface relaxes to a steady state shape, consisting of a flat height profile that couples to a curved ridge near the contact line.
Profiles ͑a͒-͑e͒ in the figure correspond to E = 0°and increasing Ca.The effect of Ca is to increase the curvature of the ridge, which in turn makes D larger.This is also observed for runs ͑f͒-͑j͒ and ͑k͒-͑o͒, which correspond to E = 25°and E = 90°, respectively.Measured contact angles are shown in Fig. 4 as a function of Ca.As seen in the plot, the profile changes from a wedge to a nose configuration by increasing the capillary number, even for E = 0°.As observed in experiments by Veretennikov et al., 15 the nose configuration has a multivalued height profile, a condition that is fulfilled only for D Ͼ 90°.For sufficiently large Ca, we achieve contact angles as large as D Ӎ 140°, a situation that renders the effective interaction between the wall and the fluid superhydrophobic.
Looking again at Fig. 3 we notice that Ca has an effect on the size of the capillary ridge, given by the maximum film thickness h max and by the width over which h is larger than h c ͑see Fig. 1͒.This effect is clearly different for hydrophilic and hydrophobic substrates.For a given Ca value, we observe thicker and wider ridges for hydrophobic substrates; runs ͑a͒, ͑f͒, and ͑k͒ were carried on at the same Ca value, still h max is approximately twice as large in run ͑k͒ than in runs ͑a͒ and ͑f͒.This tendency persists regardless of the value of the capillary number.Figure 5 shows plots of h max as a function of Ca for each contact angle considered.For E = 0°and E = 25°͓Figs.5͑a͒ and 5͑b͔͒, h max has a weak dependence with Ca, while for E = 90°͓Fig.5͑c͔͒, we find that h max decreases strongly as Ca increases.This effect is understood if we consider the limiting case of Ca→ 0, where the contact line is expected to recede due to an uncompensated Young force.As the contact line recedes the ridge grows without bound due to mass conservation. 2A finite Ca adds a driving force to the picture.The ridge then grows until the driving force and the Young force is balanced by viscous stresses.Given that the Young force is larger for larger E , the size of the ridge is expected to be larger for hydrophobic substrates.
We draw as a conclusion that the wetting properties of the substrate are relevant for the shape of the front.In contrast with lubrication theory, D is not the only control parameter of the system.One can already anticipate the relevance of this observation to the stability of the contact line.If the instability is controlled by the size of the ridge, then hydro-phobic substrates will, in general, give rise to unstable fronts, while the opposite will be expected for hydrophilic substrates.

B. Inclined substrates
To study the effect of the inclination angle ͑D 0͒ on the shape of the front we perform runs at five inclination angles for E = 90°at fixed Ca. Figure 6 shows stationary profiles for these runs.The main effect of decreasing the inclination angle is to flatten the profile.In turn, the dynamic contact angle decreases with decreasing ␣.
We extend our simulations for various Ca and E values.Figure 5 shows plots of h max / h c as a function of Ca for each ␣ and E considered.The effect of wetting properties is diminished as the inclination angle increases.For very small ␣, h max / h c → 1, regardless of the capillary number and equilibrium contact angle.Figure 7 shows plots of the excess contact angle, D − E , as a function of Ca.In all cases the contact angle increases with increasing force, as expected.The slope of the plots decreases for large Ca, indicating that the contact angle saturates to a limiting value.The effect of decreasing the inclination angle is to make the excess contact angle smaller for a given Ca value.

VI. LINEAR STABILITY ANALYSIS
Having explored the effect of wetting properties in the unperturbed front, we explore how these affect the evolution of a slightly perturbed contact line.Here, simulation parameters are chosen so the dynamic contact angle lies not only in the hydrophilic regime, but also in the hydrophobic and superhydrophobic regimes, which have not been studied previously.Front destabilization occurs by virtue of a spatial variation of the capillary ridge.When the front is perturbed transversely ͑see Fig. 1͒, the ridge at the troughs of the perturbation relaxes due to capillarity, generating transverse flows that feed the peaks.As a consequence, the peaks increase their local thickness and, in turn, their velocity.This is countered by surface tension, which tends to flatten the perturbed contact line.For long wavelengths surface tension is outweighed by the driving force and the front destabilizes.

A. Vertical substrates
We first consider a vertical substrate and explore the dynamics for five different equilibrium contact angles: E = 0°, E = 25°, E = 60°, E = 75°, and E = 90°at Ca= 0.41.In steady state, the dynamic contact angle for each substrate is D = 71°, D = 78°, D = 97°, D = 106°, and D = 128°, respectively.For a given steady front, we impose a single-mode perturbation to the contact line, x P = x U + ␦ cos͑ky͒, where x U and x P are the unperturbed and perturbed x-coordinates of the contact line, respectively, k is the perturbation mode, and ␦ the perturbation amplitude.In the linear regime, ␦ ϳ exp͑t͒, where is the growth rate of the perturbation.To estimate we measure ␦ as a function of time and perform a linear fit to the log͑␦͒-t data.
Figure 8͑a͒ shows plots of / c as a function of k / k c , where c =1/ t c and k c =1/ x c are the characteristic scales of the problem.For each angle, there exists a band of unstable modes, which broadens as one increases D .There is hence a clear dependence of the stability of the front on the wetting properties of the solid.Nevertheless, the qualitative form of the dispersion relation remains the same for all substrates.In Fig. 8͑b͒ we scale the dispersion relations, scalings being → / max and k → k / k n .Here max and k n are the maximum growth rate and the first unstable mode of each dispersion relation.The collapse is remarkably good.We draw as a conclusion that the stability of wedge-and nose-shaped fronts is determined by the same mechanism.Given that the instability is driven by gravity, thicker ridges should lead to a wider band of unstable modes and a larger growth rate.This is consistent with the fact that the size of the ridge, h max , is larger for hydrophobic substrates at fixed Ca as explained in Sec.V ͑see Fig. 3͒.A signature of this effect can be appreciated in the inset of Fig. 8͑b͒, where we plot ͑k max / k c ͒ / ͑h max / h c ͒ 2 and ͑k n / k c ͒ / ͑h max / h c ͒ 2 as a function of ͑h max / h c ͒ 2 .Using this scaling, the most unstable mode shows no dependence on h max , a condition that can only be fulfilled if ͑k max / k c ͒ϳ͑h max / h c ͒ 2 .Hence, as the size of the capillary ridge increases, the front has a smaller most unstable wavelength.
Previous results obtained by Spaid and Homsy 12 correspond to the CaӶ 1 ͑and therefore D Ӷ 1͒ limit.To regularize the contact line singularity they introduce a slip length at the solid, which affects the stability of the front.For tan͑ D ͒ Ͻ 0.1 they find a most unstable mode k max / k c = 0.5, the corresponding growth rate lies approximately in the range of 0.05ഛ max / c ഛ 0.2, and a first unstable mode in the range of 0.7ഛ k n / k c ഛ 0.9.In their simulations, larger max and smaller k n are observed as the slip length becomes small.To recover these results we fix E = 0°and Ca= 0.041.The dynamic contact angle of the steady front is D = 28°.The resulting dispersion relation is shown in Fig. 8͑a͒ and should be compared to the one corresponding to D = 71°w hich was done at Ca= 0.41 and equal E .The main effect of decreasing the capillary number is to increase the growth rate of the perturbation.This is because the size of the ridge ͑here given by its width͒ is typically larger for small Ca ͓see Figs.3͑a͒-3͑e͔͒.For Ca= 0.041 we obtain k max / k c Ӎ 0.5, k n / k c Ӎ 0.7 and a maximum growth rate max / c Ӎ 0.07.These values agree with those reported by Spaid and Homsy.Hence our results converge to the lubrication theory results at low Ca.As expected, this low-Ca dispersion relation is qualitatively equal to the rest as it collapses to the same universal curve shown in Fig. 8͑b͒.

B. Inclined substrates
We now consider the effect of the substrate inclination angle.We fix E = 90°and Ca= 0.41 and consider five inclination angles, ␣ = 60°, ␣ = 45°, ␣ = 30°, ␣ = 15°, and ␣ = 5°.As shown in Fig. 6, the profile gradually flattens as the inclination angle is decreased.For ␣ = 5°the bump disappears completely.Dispersion relations for each front are shown in Fig. 9.The effect of the normal force, in agreement with previous results, 11,13 is to suppress the instability.Indeed, for ␣ Ͻ 15°all wavelengths considered are stable.The maximum height of the front at this inclination angle is h max / h c Ӎ 1.1 The critical inclination angle for the front destabilization is expected to increase for wetting fluids.To illustrate this effect, consider, for instance, the low-Ca region of the maximum height plots shown in Figs.5͑a͒-5͑c͒ for, say, ␣ = 30°.At this inclination angle and for Ca= 0.08, we have h max / h c = 1.02, h max / h c = 1.15, and h max / h c = 1.63 for E = 0°, E = 25°, and E = 90°, respectively.Therefore, the front should be stable for E = 0°as the height approaches the limiting value h / h c → 1.For E = 25°, the ridge must be marginally unstable, given that the bump size is close ͑and above͒ the threshold value obtained for E = 90°.For E = 90°, the front must be unstable.Therefore, the critical destabilization angle increases as E decreases.From Fig. 5, the critical inclination angle to destabilize the front at E =0°is ␣ c Ӎ 60°.

A. Sawtooth and finger patterns
Having explored the linear stability of the front, we address the long time behavior of the instability for the case of a single mode.We are interested in the effect of the wetting properties on the emerging patterns.We consider four equilibrium contact angles, namely, E = 25°, E = 60°, E = 75°, and E = 90°.For each contact angle, we consider domains of a size comparable to the most unstable wavelength predicted by the linear stability analysis.

Vertical substrates
We start by considering vertical substrates.For each equilibrium contact angle we perform runs at different Ca.After the destabilization of the contact line, we observe a transient in which a single pattern emerges.Eventually, the pattern acquires a steady shape and propagates with constant velocity.
In Fig. 10 we show the evolution of contact line profiles for E = 25°and E = 90°at different capillary numbers.The evolution of the interface for angles in between is similar.Runs ͑a͒-͑c͒ in the figure correspond to E = 90°and increasing Ca.They exhibit a finger pattern that propagates at constant velocity.Wider patterns are obtained for runs ͑d͒-͑f͒, which correspond to E = 25°and increasing Ca.Of these, runs ͑d͒ and ͑e͒ evolve toward a steady finger at long times, while run ͑f͒ saturates to a steady length.This corresponds to the sawtooth shape reported in previous studies. 7,8,13able II summarizes Ca values for each of these runs.It is clear from the data that the pattern is not entirely determined by Ca, as both wide and thin patterns are observed for the same range in Ca.To characterize these patterns we measure their length L defined as the distance between the minimum ͑base͒ and maximum ͑tip͒ x-coordinates of the contact line.From Fig. 10, L is strongly affected by the wetting properties of the fluid.In Fig. 11 we show plots of L as a function of time for runs ͑a͒-͑f͒.For runs ͑a͒-͑c͒, the time variation of L is larger than for runs ͑d͒-͑f͒.In fact, for run ͑f͒, L is constant over time, meaning that the sawtooth pattern saturates to a given length in the nonlinear regime.To quantify this effect we measure the growth rate L ˙, which we summarize in Table II.For both equilibrium contact angles, the growth rate decreases as Ca increases.
Figure 12 shows three-dimensional plots of the free surface when the L ˙has reached a constant value.We observe that larger patterns are obtained when the difference between the size of the capillary ridge at the tip and the base is also large.From lubrication theory, the local velocity is expected to be u ϳ h max ͑y͒ 2 and hence L ˙͑t c / x c ͒ϳ⌬, with ⌬ = ͓h max ͑y A ͒ 2 − h max ͑y B ͒ 2 ͔ / h c 2 .Here, subscripts A and B stand for the positions of the tip and the base of the pattern ͑see Fig. 1͒, Fig. 13 shows a plot of L ˙͑t c / x c ͒ as a function of ⌬.The growth rate follows a single curve when plotted as a function of this quantity, regardless of the capillary number and equilibrium contact angle.Below a threshold value, ⌬ 0  = 0.33Ϯ 0.09, the growth rate vanishes, meaning that the amplitude of the pattern saturates while there is still a finite difference of the ridge thickness between its tip and its base.This effect is a consequence of the balance between surface tension, which tends to flatten the contact line and the uneven mass distribution caused by the instability.

Inclined substrates
We now study the effect of the inclination angle in the shape of the pattern.We consider E = 90°and Ca= 0.41 and perform simulations at ␣ = 90°, ␣ = 60°, and ␣ = 45°.Figure 14 shows the evolution of the contact line as ␣ is varied.The effect of decreasing ␣ is to decrease the growth rate and to increase the width of the pattern.The relevant effect of ␣ is to flatten the free surface, a situation that is also attained by increasing substrate wettability.Hence, the inclination angle and substrate wettability have similar effects.Indeed, the pattern growth rate in these runs falls in the same curve when plotted against ⌬, as shown in Fig. 13.

B. Finger dynamics
In the following we investigate whether the same finger solution is always observed regardless of the initial state of the system.We first examine the evolution of various initial conditions composed of a single mode.Next we consider the evolution of a multiple mode perturbation, where the most unstable mode is expected to dominate the dynamics according to linear stability results.Finally, we examine the effect of the system size.To compare results in this section to the single finger solution parameters values are fixed according to run ͑d͒ in Table II, which corresponds to a finger in the hydrophobic regime.

Effect of initial conditions
We consider a system of size W Ӎ ⌳ max , where ⌳ max is the most unstable wavelength according to the linear stability analysis.Initial conditions are fixed as ͑i͒ a single-mode perturbation with ⌳ = W, ͑ii͒ a parabolic perturbation, and ͑iii͒-͑v͒ a fingerlike perturbation with initial finger widths = 0.45, = 0.63, and = 0.80, respectively, where is measured in units of the system size.The evolution of the contact line is depicted in Fig. 15.Run ͑i͒ shows the evolution of the ⌳ = W perturbation, which grows and saturates to the finger pattern.The next four cases ͓͑ii͒-͑v͔͒ show all the same type of evolution.The interface rapidly relaxes from its initial configuration at the early times of the simulation and saturates to a finger pattern, which for all cases has the same finger width and pattern growth rate as that from the run ͑d͒ in Table II.These results indicate that the stationary solution is not sensitive to initial conditions.For a perturbation composed of more than one mode we expect that the most unstable wavelength dominates the dynamics, according to linear stability results.We consider a system of width W Ӎ 2⌳ max and an initial perturbation composed of four modes.Corresponding wavelengths are ⌳ = W, ⌳ = W / 2, ⌳ = W / 4, and ⌳ = W / 8, of which ⌳ = W / 2 is expected to grow faster than the rest.The amplitude is set to a small value ͑about two lattice spacings͒ with random deviations uniformly distributed.As shown in Fig. 16, all wavelengths except for ⌳ = W / 2 rapidly decay.The dominant mode then grows until a two-finger solution is obtained.Each of the fingers turns out to be identical to the single finger solution obtained for run ͑d͒, their length being dictated by initial conditions.

Effect of system size
We next focus on the effect of the system size in the observed pattern.We consider four system sizes, W = ⌳ max , W =3⌳ max / 2, W =2⌳ max , and W =3⌳ max .For W = ⌳ max we fix a two-finger initial condition, while for all other runs we fix a single-mode perturbation of wavelength ⌳ = W. Figure 17͑a͒ shows the evolution of the contact line for the twofinger initial condition.Even though the amplitude of both fingers is large, the system size allows only the growth of one finger.After a transient in which the shorter finger declines, the expected one-finger solution is thus recovered.
Figures 17͑b͒ and 17͑c͒ show the evolution of the contact line for W Ͼ⌳ max .Even though a single mode was present at the initial condition, the system evolves and additional modes emerge.In all cases, two-finger stationary solutions are obtained.Fingers have the same shape, width, and growth rate as the single finger solution.Nonetheless, the distance between fingers varies sensibly with increasing system size.Such a dispersion in pattern spacing is in agreement with early observations by de Bruyn, 3 who reported deviations of as much as 25% in his measurements.This was later confirmed by Kondic and Diez, 13 who pointed out that growing fingers interact depending on the initial conditions, thus generating steady patterns spaced with a certain degree of dispersion.Fingers grow whenever a sufficient amount of mass is available, i.e., when the separation between fingers lies in the band of unstable wavelengths given by linear stability.

VIII. DISCUSSION AND CONCLUSIONS
The flow of forced thin films on dry inclined substrates has been investigated for different substrate wettabilities, as well as for a wide range of inclination angles at large capillary numbers.Film dynamics explored in this work range from the hydrophilic to the superhydrophobic regimes.
Our results indicate that the stability of the contact line is influenced by the wetting properties of the solid.Wedge shaped fronts are generally more stable than nose-shaped ones.However, the stability of noses and wedges remains qualitatively the same, as all dispersion relations collapse to a single curve.The reason for this universal behavior lies in the fact that the main role of wetting properties is to fix the size of the capillary ridge, which is the relevant parameter for the stability of the front.The more hydrophobic the substrates the larger the ridge size.Then, hydrophobic substrates tend to destabilize the front.The size of the ridge can also be controlled via the inclination angle of the substrate.Accordingly, wider bands of unstable modes are observed for large inclination angles.For inclined substrates, the normal component of gravity stabilizes the film for completely wetting fluids, in agreement with previous studies. 13,14This effect is observed for partially wetting and nonwetting fluids as well.
The longtime evolution of the instability gives rise to the growth of two types of pattern; either sawtooth shaped patterns or finger shaped ones are observed.The substrate wetting properties and inclination angle have a role in the determination of the type of pattern.Fingers are observed at either hydrophobic substrates or large inclination angles, while the opposite is found for the sawtooth pattern.A direct connection with size of the capillary ridge then follows.Hydrophobic substrates and large inclination angles favor the formation of large capillary ridges.In such case, surface tension is not enough to balance the driving force given by the weight of the ridge along the contact line.Therefore, a finger develops.On the other hand, when the substrate is hydrophilic or the inclination angle is small, the size of the ridge becomes small accordingly.Surface tension then balances the driving force in the nonlinear regime, thus giving rise to a saturated sawtooth pattern.Indeed, we find that all growth rates are parametrized by the difference of the size of the ridge between the base and the tip of the pattern.The robustness of the finger solution is studied.We have found that for a given choice of Ca, ␣, and E , the same finger solution is obtained, regardless of the initial condition and system size.Fingers grow whenever a sufficient mass reservoir is available, i.e., whenever the separation between them is comparable to the most unstable wavelength predicted by linear stability.

FIG. 1 .
FIG. 1. ͑Color online͒ Thin film flowing down an inclined plane.The height of the film h is constant and equal to h c in the upstream region while it bends into a ridge close to the contact line, where it meets the solid with a dynamic contact angle D .The maximum height of the film h max ͑y͒ determines the size of the ridge.Upon a perturbation of the contact line, the ridge varies spatially.Points A and B correspond to the peak and the trough of the perturbation.Relaxation of the ridge at point B tends to feed mass to point A, therefore generating a destabilization mechanism.

FIG. 2 .
FIG. 2. Height profile comparison between lattice-Boltzmann and lubrication theory in the small Ca and D limit.Parameter values are D = 0, Ca= 0.041, and E = 0°.͑a͒ D = 0.074 and D = 28°and ͑b͒ D = 0.006 and D = 48°.Deviations from the lubrication prediction are comparable to the size of the interface, represented by the black segment located in the ridge of each profile.

FIG. 3 .
FIG. 3. Stationary xz-interface profiles for D = 0 at various E and Ca.The interface is defined by the isosurface =0.The x-and z-scales are the same for all plots.Each column corresponds to a different equilibrium contact angle; from left to right: E = 0°, E = 25°, and E = 90°.For each column Ca increases downward ͑rows do not correspond to the same Ca value͒.

FIG. 5 .
FIG. 5. Maximum film thickness as a function of Ca. ͑a͒ E = 0°.͑b͒ E = 25°.The results are presented explicitly for ␣ and Ca, as h does not scale with the lubrication parameter D. ͑c͒ E = 90°.

FIG. 12 .
FIG. 12. ͑Color online͒ Interface profiles at different contact angles and capillary numbers.Parameters for each run are summarized in TableII.

FIG. 16 .
FIG. 16.Contact line evolution for a perturbation composed of four modes.

TABLE I .
Simulation parameters for runs presented in Fig.3.
FIG. 4. Dynamic contact angle as a function of Ca.

TABLE II .
Simulation parameters and pattern growth rate for vertical substrates.
FIG. 11.Pattern length as a function of time.Letter keys correspond to parameters shown in Table II.͑a͒-͑c͒ E = 90°and increasing Ca; ͑d͒-͑f͒ D = 25°and increasing Ca.