The impact of microscopic liquid drops on solids with a variety of surface characteristics is studied using numerical simulations. The focus is on relatively low impact velocities leading to bouncing or spreading drops, and the effects of wettability. Molecular dynamics and lattice Boltzmann simulation methods are used for nanometer-sized and continuum drops, respectively, and the results of the two methods are compared in terms of scaled variables. We consider surfaces which are flat, curved or pillared, with either homogeneous interactions or cross-shaped patterns of wettability. In most situations we observe similar drop behavior at both length scales; the two methods agree best at low impact velocities on wettable surfaces while discrepancies are most pronounced for strongly hydrophobic surfaces and for higher velocities.

The impact of liquid drops on solid surfaces is an everyday phenomenon with important applications in material transport processes in industry, and whose understanding poses challenging questions in fluid dynamics.1–3 The physical properties of the solid obviously play an important role in the behavior of impacting drops, and recent developments in fabrication techniques which provide well-characterized surface structures with specified properties allows one to study impact dynamics with much greater precision than in the past. The surfaces available for study range from chemically homogeneous to randomly disordered, from hydrophilic to super-hydrophobic, and from atomically smooth to geometrically or randomly rough.4–7 Furthermore, there are numerous possibilities for the outcome of an impact: the drop may bounce, stick, spread as a disk, spread into fingers or form a splash.8–10,12 A spreading liquid lamella may rupture during impact into satellite drops or fragments.10,13,14 If the substrate temperature is above the liquid boiling point, the drop may levitate (the Leidenfrost effect15,16) or even leave the surface entirely, depending on its size.17 A recent systematic review of the phenomena can be found in Ref. 8.

This paper presents numerical simulations of drop impact, with two principal motivations. First, we aim to predict in a general way how the behavior of an impacting drop is affected by the chemical and structural characteristics of the surface. Although slightly simplified physical models are used for surface properties, we hope to elucidate the trends and compare the results to experiments and earlier calculations where possible. We employ two distinct numerical simulation methods: molecular dynamics (MD), which accurately captures the behavior at the scale of individual atoms in modestly-sized systems, and the lattice Boltzmann method (LB), which solves the continuum equations of motion using a particle-based algorithm. Ostensibly the two methods operate at very different length and time scales, but in principle address the same phenomena. Our second goal in this paper is to compare the results of the two types of calculation in the non-trivial flow configurations arising in drop impact. An earlier comparison of the two methods for the confined (Poiseuille) flow of a wetting liquid18 found good agreement, but drop impact is a more stringent test.

To anticipate the results somewhat, we find that at low values of the Reynolds and Weber numbers (Re ≲ 20 and We ≲ 100) the two methods produce very similar drop behavior on wettable surfaces. On nonwettable surfaces and at higher impact velocities, and in particular in the splashing regime, deviations appear. In high velocity impacts, there is first the general issue that any splash produces a highly ramified fluid body and secondary droplets, and resolving this structure in terms of individual molecules requires a very large scale simulation. A further difficulty is that a drop several nanometers in diameter with high Re and We has a high velocity O(100 m/s), comparable to the atomic thermal velocity. One consequence is that the Mach number becomes significant and compressibility effects may appear, and a second is that impact produces temperatures exceeding the liquid/vapor coexistence value,10 and the liquid drop tends to disintegrate rather than splash in the conventional sense. Neither of these features is present in most experiments and continuum calculations, and we therefore focus on slower impacts.

A number of earlier papers have studied drop impact on heterogeneous surfaces of various kinds, using both experimental and numerical methods. Cross, stripe, and chessboard chemically-patterned surfaces have been studied numerically in the lubrication approximation by Schwartz and co-workers,19,20 using a geometry identical to one considered below. Experimentally, Vaikuntanathan et al.21 considered drop impact at the junction of hydrophobic and hydrophilic half-plane surfaces. Cheng22 measured the dampened vibratory motion of water drops impacting onto a coal surface both theoretically and experimentally, and detailed experimental studies of the fingering were conducted by Marmanis and Thoroddsen.12 They counted the number of fingers seen on the splatter left behind after a colored drop impact on a thick rough linen paper. Mock et al.23 studied impacts onto hydrophilic spots and arrays of spots. A pattern made of radiating thin spokes is the subject of experiments by Lee et al.24 and annular patterned surfaces are studied by Kim et al.,25 and experiments by Katsuragi26 involve a drop impact on a granular layer surface. The dynamics of water drop impact at high impinging velocity onto super-hydrophobic substrates is experimentally investigated by Tsai et al.,27 and multiscale pillar surfaces are studied by Lohse et al.6,28 and Ruckenstein et al.29,30

The paper is organized as follows. In Sec. II we review the implementation of the MD and LB techniques used in our simulations. Some of the configurations we studied revisit previous work while others are novel, to our knowledge. In Sec. III, for example, we begin with the simplest case of homogeneous surfaces of constant wettability and compare MD and LB results in some detail. Next, in Sec. IV we consider “cross patterns” involving distinct wettable and nonwettable regions organized in the shape of a plus sign. We discuss how different combinations of the wettabilities produce different morphologies in drop impact. In Sec. V, we turn to more novel rough surfaces, consisting of a periodic array of square pillars on a plane. Here we observe a variety of possible drop behaviors, due to the interplay of Cassie and Wenzel states. Finally, in Sec. VI, we address some of the effects of larger-scale roughness by considering drop impact on curved but otherwise smooth surfaces. Here we discuss the differences between concave and convex surfaces, and the effect of curvature on drop spreading. Conclusions are given in Sec. VII.

The formulation of the MD calculations in this paper is almost identical to that of a previous paper on drop impact on an homogeneous non-wetting plane,10 except that more general solid surfaces are considered.

In a typical molecular simulation, shown in Fig. 1, a drop plus vapor system of 115 392 dimer molecules floats above a flat surface composed of 245 650 wall atoms arranged into one layer of fcc cells, all placed in a three-dimensional box of size (300σ, 120σ, 300σ). The length unit σ is the core size (approximate molecular diameter) of the Lennard-Jones potential acting between all nearby atomic pairs,

\begin{equation}U(r)=4 \epsilon \left[ \left( \frac{\sigma }{r} \right)^{12} - c \left( \frac{\sigma }{r} \right)^6 \right] \qquad r \le 2.5\sigma .\end{equation}
U(r)=4εσr12cσr6r2.5σ.
(1)

Here ε is the depth of the potential well and the parameter c controls the strength of attraction between two atoms. We choose the standard value c = 1 for fluid and wall atom pairs but adjust c to vary the wettability for the fluid-solid interactions. The correspondence between the contact angle θ and c can be either roughly estimated by cos θ ≃ 1 − 2cρsl,31 where ρs, l are the liquid and solid densities, respectively, or more precisely measured by a separate simulation, with the results given in Table I. The calculations employ our own serial code, and require about a day for each impact simulation.

FIG. 1.

Snapshot of a drop and vapor above a solid surface before impact in a MD simulation.

FIG. 1.

Snapshot of a drop and vapor above a solid surface before impact in a MD simulation.

Close modal
Table I.

Contact angle θ as a function of the wettability c for a dimer system of density ρ = 0.8mσ−3 at temperature T = 0.8ε/kB.

c 0.4 0.67 0.8 0.9 1.0 
θ 180° 134° 96° 73° 50° 0° 
c 0.4 0.67 0.8 0.9 1.0 
θ 180° 134° 96° 73° 50° 0° 

Dimer molecules are implemented by linking two fluid atoms with a FENE potential,32 and periodic boundary conditions are applied in the x and z dimensions parallel to the surface. A Nosé-Hoover thermostat is used to equilibrate the system before impact at temperature T = 0.8ε/kB, but afterwards only the solid temperature is controlled while that of the fluid is allowed to vary. To initiate drop impact, the velocity of each drop molecule is given a downward shift of magnitude v0, ranging from 0.5 to 2.0σ/τ.

To quantitatively describe the drop dynamics after impact some numerical characterization of the drop shape is needed. A typical drop shape slightly after impact is shown in Fig. 2, along with its boundary. Provided that the solid surface is uniform and we avoid the splashing regime, the drop is observed to spread with approximate radial symmetry in the plane of the surface. Using a cylindrical coordinate system whose (y-)axis is a vertical line through the center of the drop, three obvious geometrical parameters characterizing the r-y profile are the drop height h, the spreading radius Rm, and the radius of the contact area Rc. The intrinsic length scales of the MD and LB calculations are different, and to facilitate a comparison we normalize all lengths by the initial drop radius R0.

FIG. 2.

Snapshot of an early stage of drop spreading (left) and the corresponding interfacial contour (right), as obtained by the method described in the text. The definitions of h, Rm, and Rc are illustrated and the dashed line is the wall.

FIG. 2.

Snapshot of an early stage of drop spreading (left) and the corresponding interfacial contour (right), as obtained by the method described in the text. The definitions of h, Rm, and Rc are illustrated and the dashed line is the wall.

Close modal

In drawing Fig. 2 the boundary of the drop is assumed to be a line, but in fact in both computations and in reality a liquid vapor interface is a transition region of finite thickness. Our procedure is based on the fluid density field ρ(r, y) in the cylindrical coordinate system above. Before impact we observe that away from the interface the density has constant values ρl and ρv in the center of the drop and in the distant vapor regions, respectively, and we identify the interface as the midpoint curve on which ρ(r, y) = (ρl + ρv)/2. Typical MD values are ρl = 0.8 mσ−3 and ρv = 10−3ρl, with an interfacial thickness O(2-3 σ). Strictly speaking, in MD simulations at any single time ρ is a sum of delta-functions centered at the current atomic positions but if we time-average over a 1 τ interval a smoothly-varying field results. In the LB the situation is simpler because this continuum density field is a direct output of the calculation.

While h and Rm are now well defined, there is some further ambiguity in the measurement of Rc because “contact” is ill-defined at molecular scales. There is always a finite distance between the adjacent liquid and wall atoms (due to the repulsive core of the potential), whose value depends on the wettability of the surface. Our rule is to choose a critical value bc of the gap between the interface contour defined above and the average position of the top wall atoms, and say that contact occurs where the gap b(r) < bc.

In analyzing the time dependence of the drop parameters a starting time must be chosen, which is also ambiguous because of the finite range of the interaction. Our convention is to start the measurement at the moment when the drop “geometrically” touches the surface, meaning when when h = 2R0. In the MD simulations, at this time there are always drop molecules within the critical gap thickness and contact in the sense of the previous paragraph has been made, but in the LB case we find that the drop height h is distorted by up to 5% before the contact line forms. In experiment an additional complication may arise–an air bubble trapped in the center immediately after drop impacts the wall,33 which may delay the formation of a contact line, but such bubbles are not observed here.

The LBM was first used by Dupuis and Yeomans34 to model drops on a super-hydrophobic surface. This technique was consequently used to study Cassie-Wenzel transition,28,35,36 and drop motion on super-hydrophobic surfaces.37–39 Recently, Connington and Lee40 presented an extension of the method of Lee and Liu41,42 to implement boundary conditions for the complex geometry of super-hydrophobic surfaces. The model is validated by comparing with experimental work of Barbieri et al.43 It was found that LB simulation reproduces the experimental contact angles well, both in the Cassie and in the Wenzel state, but the transition between states is not accurately captured, which is due to the diffuse interface detecting an initial contact with the substrate before a sharp interface would, initiating the transition precipitately. When the sharp interface limit is approached, it was found that the transition is captured more accurately.

In the lattice Boltzmann method, one follows the evolution of a density distribution function for fictitious particles moving on a lattice, which depends on position and velocity. The velocity is discretized so as to have one lattice spacing per time step, and particle motion consists of one streaming step between neighboring lattice sites followed by a collision step, treated in the relaxation-time approximation. In order to model two-phase flows intermolecular interaction forces are incorporated into the governing equations for the particle distribution functions. In this study we employ the scheme developed by Lee and Liu41,42 using a scalar order parameter C, which satisfies the convective Cahn-Hilliard equation ∂tC + ∇ · (uC) = M2μ, involving an adjustable mobility M and a chemical potential μ. The latter is given by the derivative of the free energy with respect to the order parameter

\begin{equation}\Psi =\int _S\left[ E_0(C)+\frac{\kappa }{2}|\nabla C|^{2}\right]dV +\int _S\left[\phi _0-\phi _1C_s+\phi _2C_{s}^{2}-\phi _3C_{s}^{3}+ \ldots \right]dS,\end{equation}
Ψ=SE0(C)+κ2|C|2dV+Sϕ0ϕ1Cs+ϕ2Cs2ϕ3Cs3+...dS,
(2)

involving a bulk free energy E0(C) = βC2(C − 1)2, where β is a constant, a gradient term whose coefficient κ controls the surface tension and interfacial thickness, and surface terms which control the solid-liquid interactions, involving the surface concentration Cs. Proper control of the wetting properties of a fluid on a solid substrate is obtained by formulating the wall free energy polynomial boundary conditions (linear, quadratic, and cubic).44,45 In the current paper the integral of the free energy on solid boundaries employs a cubic boundary condition, where the interactions between the solid and bulk fluids are neglected and only the interactions between the solid and interface are considered. This boundary condition is appropriate for our simulations, when there is a large difference in density between the two fluids. This assumption specifies the parameter set: ϕ0 = ϕ1 = 0, ϕ2 = 1/2ϕc, and ϕ3 = 1/3ϕc, where ϕc is a constant to be chosen to recover the desired contact angle at equilibrium. The dimensionless wetting potential

$\Omega _c=\phi _c \sqrt{2\kappa \beta }$
Ωc=ϕc2κβ⁠, is related to the equilibrium contact angle by cosθeq = (σsg − σsl)/σ = −Ωc. It was concluded that the equilibrium contour is preserved very well on a wetting surface, but there is a noticeable reduction in the radius of the drop on a non-wetting surface. The proposed boundary conditions were shown to effectively eliminate spurious currents at equilibrium, and the total mass in the system does not change by more than 1% of its initial value. They all are also able to produce almost analytical solution of the equilibrium contact angle and density distribution of a static drop on solid surfaces with different wetting characteristics.44 The liquid/vapor interface profile is therefore found by minimization of the free energy, with cubic boundary conditions at a solid surface.

To illustrate a typical LB calculation, we consider a micron-scale drop impact with liquid/vapor density ratio of 842 and viscosity ratio of 51, and with different wettabilities on flat dry surfaces.41 We refer to micron-scale here because the method operates at continuum length scales, greater than a nanometer at least, and gravitational forces are omitted, as would be appropriate in this regime. A three dimensional liquid water drop on a solid surface with equilibrium contact angles ranging from 0 to 180 is generated at the bottom center of a 135 × 135 × 112 computational domain for a D3Q27 lattice. The boundaries are all symmetric except at the solid surface, where the wall boundary condition44 is imposed. The interface thickness and initial drop radius are ζ = 5 and R0 = 37.5, respectively, where all parameters are given in lattice units. The impact process is described by three independent dimensionless numbers, the contact angle, Weber number (We), and either the Ohnesorge number (Oh) or the Reynolds number (Re), defined as follows:

\begin{equation}We=\frac{\rho _1{u_0}^{2}D_0}{\sigma } \quad \mbox{\rm and}\quad Oh=\frac{\eta _1}{\sqrt{\rho _1\sigma D_0}} \quad {\rm or}\quad Re=\frac{\sqrt{We}}{Oh} ,\end{equation}
We=ρ1u02D0σandOh=η1ρ1σD0 or Re=WeOh,
(3)

where u0 is the drop impact speed, D0 is the drop diameter prior to impact, and η1 and ρ are the fluid's viscosity and density. (Note that in the previous paper,10 Re and We are defined differently, using the drop radius instead of the diameter.) We consider Oh = 0.477 and the two Weber numbers, 128 and 288, for a fixed density ratio of 513 (i.e., η1 = 1.0 and η2 = 1.95 × 10−3), and a viscosity ratio of 14. The grid dependency of the results is tested by repeating the calculations using systematically refined 2-dimensional grids for a D2Q9 lattice, with three different grid resolutions: 90 × 90, 135 × 135, and 202 × 202, corresponding to droplet radii of 25, 37.5, and 56.25, respectively. The time evolution of the dimensionless diameter or spreading ratio ξ* = D/D0 for these three grid sizes is shown in Figure 3 for We = 128, Re = 24 with equilibrium contact angle of 78°. From the figure we conclude that the results are stable once the grid size reaches 135 × 135, and this resolution is used in the subsequent calculations.

FIG. 3.

LB validation. (Left) Mesh dependency results: spreading ratio vs. time for the three resolutions indicated. (Right) Maximum diameter of the spreading drop, normalized by the drop radius, as a function of the Weber number. The solid line indicates slope 1/4.

FIG. 3.

LB validation. (Left) Mesh dependency results: spreading ratio vs. time for the three resolutions indicated. (Right) Maximum diameter of the spreading drop, normalized by the drop radius, as a function of the Weber number. The solid line indicates slope 1/4.

Close modal

To validate the computations we compare the results with experimental data by Clanet et al.47 for water and mercury drop impact on a super-hydrophobic surface with equilibrium contact angle of θeq = 180°, a fixed density ratio of 842, and a viscosity ratio of 51. The Ohnesorge number is 0.015, and the Weber number varies between 8 and 90. The data are consistent with a scaling law for the maximum spreading ratio, DmaxD0We1/4, which is as a solid line shown in Figure 3, along with the corresponding LB results. The agreement between experiment and LB calculation in this situation is quite good. The prefactor coefficient for the solid line equals to 0.9, while this value for LB simulation and water and mercury droplets is 1.1, 0.887, and 0.89, respectively.

We begin with the simplest case—impact on a flat and homogeneous solid surface. A direct comparison between the (azimuthally-averaged) contours of spreading drops at Re = 24 and We = 128 at two contact angles is given in Fig. 4. The drop shapes are very similar at the earlier times after impact, for times below t* = R0/v0, and begin to differ afterwards although there are no gross differences. The shape very near the contact line, however, is somewhat different throughout the process, which presumably reflects differences in the modeling of the solid/liquid interactions.

FIG. 4.

Azimuthally-averaged drop shapes during impact for Re = 24 and We = 128. Top: contact angle θ = 96°, bottom: θ = 50°; MD is solid line (blue), LB is dashed line (red). From left to right: t/t* = 1, 2, 4, 8. The length unit is the initial drop radius R0 and the time unit is t* = R0/u0.

FIG. 4.

Azimuthally-averaged drop shapes during impact for Re = 24 and We = 128. Top: contact angle θ = 96°, bottom: θ = 50°; MD is solid line (blue), LB is dashed line (red). From left to right: t/t* = 1, 2, 4, 8. The length unit is the initial drop radius R0 and the time unit is t* = R0/u0.

Close modal

For a more quantitative comparison, we consider the variation of drop height and radius with time. First, in Fig. 5 we plot the height of the center of the drop in the two methods for four impact velocities at a 50° contact angle. Lower angles behave somewhat similarly, but there are substantial discrepancies in strongly non-wetting case. In the figures, the drop height h is scaled by the initial drop radius R0 and time is scaled (inertially) by the free-space transit time across a radius, t* = R0/v0. Generally, the height decreases linearly until impact, then less rapidly as the drop is squeezed against the solid, and then a weak nearly-linear rise as the drop re-expands upwards to its final state. The two methods are in reasonable agreement at the lower impact speeds but do not agree at higher Re. Next, we compare the results for different contact angles at a fixed impact velocity Re = 24 and We = 128. We see in Fig. 6 that as a function of contact angle, the h(t) curves for the two methods are similar for low contact angles but deviate as θ increases and the wettability decreases. For a completely non-wetting surface (θ = 180°) where the disagreement is largest, the MD drops simulated here have no attractive interactions at all with the solid surface and tend to float off unless there is a significant impact velocity, whereas this property is absent in the LB case. Note that the scaling behavior in these results is rather erratic: the curves for different impact velocities overlap only at early times, while those for different contact angle in wetting and partially-wetting cases form a nearly-universal curve in the LB case but do so only approximately in MD.

FIG. 5.

Time dependence of the drop height after impact for various impact speeds at contact angle 50°. Left: MD; right: LB. ○: Re = 10; △: Re = 19; □: Re = 23.7; ▽: Re = 35.6.

FIG. 5.

Time dependence of the drop height after impact for various impact speeds at contact angle 50°. Left: MD; right: LB. ○: Re = 10; △: Re = 19; □: Re = 23.7; ▽: Re = 35.6.

Close modal
FIG. 6.

Time dependence of the drop height after impact for various wettabilities. at Re = 24. Left: MD; right: LB. Contact angle: ○: θ = 180°; □: θ = 96°; △: θ = 50°; +: θ = 0°.

FIG. 6.

Time dependence of the drop height after impact for various wettabilities. at Re = 24. Left: MD; right: LB. Contact angle: ○: θ = 180°; □: θ = 96°; △: θ = 50°; +: θ = 0°.

Close modal

The variation of drop radius with time is consistent with this behavior. Since the drop shapes are irregular and not simple functions there is no unique characterization of their radius, so we have examined the two obvious choices—the maximum value Rm, corresponding to the drop size as projected from above, and the contact radius of liquid on the surface Rc, which might be viewed from below in suitable experiments on transparent solids. The two radii are most sensitive to different physical effects. The maximum radius begins at the initial value R0 at contact, increases due to inertia, slows due to viscosity, and may withdraw due to surface tension, although this depends on the surface properties. The contact radius is very sensitive to the liquid/solid interaction, and tends to be less than the maximum radius on nonwetting surfaces. It is difficult to anticipate the relative values of radii on partially-wetting surfaces and crossovers may occur. The dynamic contact angle itself is generally quite different from the equilibrium angle, since it starts at 180° at impact and tends to relax to the equilibrium angle when spreading halts.

We begin with Fig. 7 which shows the variation of the two radii with velocity for a weakly-wetting surface at contact angle 50°. In the MD case drops at all impact velocities spread and then contract, consistent with approximate incompressibility and the previous observation that the drop height rises after impact. In microscopic terms, the contraction results from weak solid/liquid attraction being overcome by liquid/vapor surface tension, whereas in continuum terms the inertially-driven spreading (Fig. 4) drives the contact angle to a value above the equilibrium angle and the drop withdraws to return to the proper value. In contrast, on a completely wettable surface, both radii increase indefinitely as the drop attempts to cover the surface, although rather slowly after initial impact on the time scale of these simulations. In general, a receding stage may or may not be observed, depending on the surface wettability. In the LB case, a noticeable withdrawal of the drop after initial spreading is seen only in the maximum radius, and there the effect is rather weaker than in MD. The contact radius hardly withdraws at all, suggesting differences in the LB and MD treatment of surface interactions. Furthermore, as in the height variation, there is a rough agreement between the two methods at least at early times only at lower impact velocities.

FIG. 7.

Time dependence of the drop radii after impact, for fixed equilibrium contact angle (50°) and various impact speeds. (a) Rm in MD; (b) Rm in LB; (c) Rc in MD; (d) Rc in LB. ○: Re = 10, We = 23; △: Re = 19, We = 82; □: Re = 23.7, We = 128; ▽: Re = 35.6, We = 288.

FIG. 7.

Time dependence of the drop radii after impact, for fixed equilibrium contact angle (50°) and various impact speeds. (a) Rm in MD; (b) Rm in LB; (c) Rc in MD; (d) Rc in LB. ○: Re = 10, We = 23; △: Re = 19, We = 82; □: Re = 23.7, We = 128; ▽: Re = 35.6, We = 288.

Close modal

A possible source of discrepancy at higher velocities is in the treatment of slip. In MD, the contact line motion is that of the atoms in that region, as controlled by the interplay of interactions with the solid and with the rest of the liquid. It is well known that the resulting slip length increases systematically as the fluid velocity increases and as wettability decreases.11 In the highest-velocity impact (Re = 35.6) substantial slip is observed on non-wetting surfaces10 which promotes rapid spreading. In the LB method, the surface interactions are constructed so as to produce a non-slip boundary condition for the velocity, but in fact the motion of the contact line is due to diffusion of the Cahn-Hilliard concentration field. (It has been shown48 that reducing interfacial thickness while keeping the mobility constant gives the sharp interface limit with a diffusion length scale corresponding to the slip length, and in agreement with the asymptotic analysis of Cox49). The resulting value of the contact line velocity differs in the two methods, as illustrated in Fig. 8. In both models, the contact line velocity is obtained from the time-dependence of the concentration midpoint value in the row of bins (sampling and computational in the MD and LB methods, respectively). Except at early times, the contact line velocity is somewhat larger in MD, which enhances the spreading. Note also the anomalous behavior in MD in the nonwetting 180° case, which shows distinctly more slip than for the other contact angles.

FIG. 8.

Time-dependence of the contact line velocity for unit impact velocity, comparing MD and LB results for different contact angles as indicated in the inset at Re = 35.6, and We = 288.

FIG. 8.

Time-dependence of the contact line velocity for unit impact velocity, comparing MD and LB results for different contact angles as indicated in the inset at Re = 35.6, and We = 288.

Close modal

There are also distinctions in the drop shapes between the two methods, arising from the fact that the MD simulations operate near the boundary of the splashing regime whereas the LB calculations produce smooth and stable drop shapes in these conditions. At the higher Reynolds numbers, the MD drops exhibit features which might be thought of a precursors to splashing: an irregular edge and transient holes appear in the spreading lamella during the expansion stage, as shown in Fig. 9. Furthermore, the drop at Re = 36 on a nonwetting surface does in fact splash, and the results for this case in the figures above are restricted to times before the splash.

FIG. 9.

Top view of the MD drop shape at Re = 36 and We = 228 for different wettabilities at t = 6t*. Left to right: θ = 0°, θ = 95°, and θ = 180°.

FIG. 9.

Top view of the MD drop shape at Re = 36 and We = 228 for different wettabilities at t = 6t*. Left to right: θ = 0°, θ = 95°, and θ = 180°.

Close modal

A second source of the discrepancy between the MD and LB calculations as impact velocity increases is the effect of Mach number. At Re = 10, in these MD simulations the Mach number Ma ≈ 0.1,10 which at first glance is worrisome for comparisons to normal liquids in typical impact experiments and simulations, but from the average density output of the calculation (the volume enclosed by the interface contours) we estimate that the maximum variation in drop volume is at most only about 2%—see Fig. 10. However, the Mach number increases linearly with impact velocity and correspondingly the volume variation rises to about 7% when Ma = 0.24, which is a cause for concern. In the LB calculations, the Mach number may be varied to some degree to address this point, but beyond Ma ≈ 0.1 the calculation fails to converge properly. However, if we extrapolate the LB calculations to the values of Ma relevant to the MD simulations, then as shown in Fig. 10 the disagreement in maximum spreading radius is significantly reduced.

FIG. 10.

Mach number effects in the LB calculations. (Left) Volume change during impact. (Right) Variation of the maximum spreading factor with Ma. The curves are computed for Ma ⩽ 1 and linearly extrapolated to larger values.

FIG. 10.

Mach number effects in the LB calculations. (Left) Volume change during impact. (Right) Variation of the maximum spreading factor with Ma. The curves are computed for Ma ⩽ 1 and linearly extrapolated to larger values.

Close modal

The wettability dependency of the spreading radii is shown in Fig. 11. Note the anomalous behavior of the fully non-wetting θ = 180° case: except at the earliest times, the variation of the radii is entirely different from the other cases where there is some degree of wetting. Furthermore, the contact radius varies erratically in the MD plot, while the maximum radius grows very rapidly and then decays very rapidly, corresponding to a drop which is barely in contact with the solid surface while rapidly expanding and then contracting from a lamella state. The LB results partly echo this behavior, although contact with the solid is maintained so that Rc varies smoothly and time-dependence of Rm is less rapid. In the other cases, there is approximate agreement between the two methods and some uniform scaling behavior at early times up to 2-3t* suggesting that inertia dominates the spreading dynamics there. Notice that the radii grow at long times in the two wetting cases (θ = 0° and 50°) but decay after initial growth in the non-wetting case (θ = 90° and 180°). Furthermore, the scaled time at which the two radii reach their maximum values in non-wetting situations is in the range 1.8–3.9 in Figs. 7 and 11, consistent with the data of Dong et al.1 (note that a different definition of scaled time is used in this paper). At early times the contact radii are independent of wettability, except for the 180° case. Similar results were found in a MD study of cylindrical drops placed at rest near a solid surface by Winkels et al.,46 who found a t1/2 growth law for the contact radius. Our results are consistent with this behavior.

FIG. 11.

Variation of spreading radii with wettability at Re = 24 and We = 128. (a) Maximum radius from MD; (b) maximum radius from LB; (c) contact radius from MD; (d) contact radius from LB. ×: θ = 180°; +: θ = 96°; △: θ = 50°; □: θ = 0°.

FIG. 11.

Variation of spreading radii with wettability at Re = 24 and We = 128. (a) Maximum radius from MD; (b) maximum radius from LB; (c) contact radius from MD; (d) contact radius from LB. ×: θ = 180°; +: θ = 96°; △: θ = 50°; □: θ = 0°.

Close modal

Finally, we compare the present simulations to results in the literature. There are direct experimental results on hydrophobic surfaces by Clanet et al.,47 which indicate that the We number dependency of the maximum spreading factor is ξmaxWe1/4,47 where ξmax is the maximum spreading diameter during impact scaled by the initial drop diameter D0. An empirical extrapolation from experiments by Scheller and Bousfield,50 a combination of data and theory by Roisman51 and models based on energy balance arguments by Pasandideh-Fard et al.52 extrapolated to the fully hydrophobic case show slightly different behavior, a power law with a different prefactor and a smaller exponent. Recent experiments by Visser et al.53 are with the model of consistent with the latter model. Our MD results in Fig. 12 lie in between the data, and exhibit qualitatively consistent behavior: we find almost the same exponent (0.24) at lower We, where the calculation is in the best correspondence to experiment (lower Ma) but a more rapid variation with exponent around 0.5 at higher impact velocities: For We > 120, ξmax scales as We0.49 for Dm and We0.56 for Dc (the two dashed lines in the figure). At still higher We the spreading lamella in the MD simulations begins to break up at its edges (see Fig. 9), and the results there are uncertain. The LB results instead vary as a single power-law which is bracketed by the experimental points and overlaps the MD points at low We. For the weakly wetting case θ = 96° there are only extrapolated or theoretical previous results, which again have some spread. The MD simulations again show distinct power-laws at low and high We, and are only in rough agreement with previous work at the higher We values. The LB points again have uniform power-law behavior and have the same relation to the other data as in the nonwetting case.

FIG. 12.

Weber number dependency of the maximum spreading factor, for contact angles left: 180° and right: 96°. MD results: ▲: ξmax calculated from Dm; ▼: ξmax calculated from Dc. LB results: △: ξmax calculated from Dm; ▽: ξmax calculated from Dc. Experiments and extrapolations thereof: |$\; \smash{\lower0.2pt\hbox{- }}\hspace*{-5.3pt}\ast\; $|∗̶: Clanet et al.;47|$ \; \smash{\lower0.2pt\hbox{- }}\hspace*{-6.4pt}\times $|☓̶: Scheller and Bousfield;50|$ \; \smash{\lower0.2pt\hbox{\hbox{---} }}\hspace*{-7pt}+\; $|: Roisman;51|$ \; \smash{\lower0.2pt\hbox{\hbox{---} }}\hspace*{-7pt}\square\; $|⬜̶: Pasandideh-Fard et al.;52|$ \; \smash{\lower0.2pt\hbox{- }}\hspace*{-5.8pt}\lozenge\; $|◊̶: Chandra and Avedisian.54 

FIG. 12.

Weber number dependency of the maximum spreading factor, for contact angles left: 180° and right: 96°. MD results: ▲: ξmax calculated from Dm; ▼: ξmax calculated from Dc. LB results: △: ξmax calculated from Dm; ▽: ξmax calculated from Dc. Experiments and extrapolations thereof: |$\; \smash{\lower0.2pt\hbox{- }}\hspace*{-5.3pt}\ast\; $|∗̶: Clanet et al.;47|$ \; \smash{\lower0.2pt\hbox{- }}\hspace*{-6.4pt}\times $|☓̶: Scheller and Bousfield;50|$ \; \smash{\lower0.2pt\hbox{\hbox{---} }}\hspace*{-7pt}+\; $|: Roisman;51|$ \; \smash{\lower0.2pt\hbox{\hbox{---} }}\hspace*{-7pt}\square\; $|⬜̶: Pasandideh-Fard et al.;52|$ \; \smash{\lower0.2pt\hbox{- }}\hspace*{-5.8pt}\lozenge\; $|◊̶: Chandra and Avedisian.54 

Close modal

As a prototypical example of a heterogeneous surface we consider drop impact on a surface with a cross-shaped region of one wettability superposed on a background with a different value. Since different regions of the edge of a drop on a cross pattern have different contact angles, its motion will no longer be radially symmetric and interesting spreading patterns result. The surface is indicated schematically in Fig. 13, where two complimentary cases are shown: a wettable cross pattern on a nonwetting substrate (case A), and a nonwetting cross pattern with a wettable background (case B). The width of the cross b is set to R0/2, and the drop impact is centered at the middle of the cross, so that the drop necessarily covers regions of both wettabilities. The equilibrium contact angles inside and outside the cross are denoted as θi and θo, respectively, and the specific values used are 73° and 180° for wetting and nonwetting regions, respectively. In MD simulations, we treat the wall atoms in different regions as different species and set the liquid-wall interaction strength c so as to produce the appropriate contact angle. In LB simulations, the contact angles are direct inputs and we simply assign different θ's to the different wall region.

FIG. 13.

Illustration of the cross patterns of wettability. Left: case A; right: case B. Dark and light regions are wettable and nonwetting, respectively.

FIG. 13.

Illustration of the cross patterns of wettability. Left: case A; right: case B. Dark and light regions are wettable and nonwetting, respectively.

Close modal

A similar problem was previously studied numerically by Schwartz et al.,19,20 although the configuration included a thin wetting layer covering the substrate to avoid a no-slip singularity in the lubrication approximation calculation. Various drop motions resulted, including spreading, migration and drop breakup, but the differences in assumptions preclude a detailed comparison with the present work.

Snapshots of the two simulations for Re = 24 and We = 128 for the two patterns are given in Fig. 14. In the MD simulations, during the spreading stage after impact, for t ⩽ 2t* in the figure, the drop spreads more slowly on wetting than on nonwetting surfaces, resulting in a rectangular shape in case A and a diamond shape in case B. When the drop recedes, it likewise does so more slowly on the wetting region, leaving a squarish shape with nodes along the wetting stripes in case A and an approximate disk with indentations along the cross in case B. The corresponding LB simulation results show a much smoother drop shape at all times, as well as a much weaker variation of the drop radius with time, but the shape are qualitatively similar in terms of protuberances and indentations along the cross.

FIG. 14.

Snapshots of the top view of the drop shape during impact at Re = 24 and We = 128. First row: MD for case A; second row: LB for case A; third row: MD for case B; fourth row LB for case B. Left to right: t/t* = 2, 4, 6, 10.

FIG. 14.

Snapshots of the top view of the drop shape during impact at Re = 24 and We = 128. First row: MD for case A; second row: LB for case A; third row: MD for case B; fourth row LB for case B. Left to right: t/t* = 2, 4, 6, 10.

Close modal

To quantify the shape for this pattern geometry, we measure the spreading radii of the drop along the four axes and four diagonal (45°) directions. Averaging along symmetry directions yields two radii, one along the axes (the cross) and the other along the diagonals (outside the cross), and their time evolution is compared to that on homogeneous surfaces with the same wettabilities in Fig. 15. In both cases, the early stages of spreading are identical to that on a homogeneous wall, but differences occur when the drops recede. In case A in the MD simulations, the drop recedes more rapidly along the wettable cross then for a homogeneously wettable surface, because it is pulled inwards by the rapidly-receding fluid in the nonwettable diagonal regions. In the corresponding LB figures the radii hardly differ until times around 3t*, and then the behavior on the two regions is close to their respective homogeneous cases. In case B, the MD fluid along the nonwetting cross is slightly slowed by the slow-moving fluid in the wetting diagonals, and the effect is weaker because the amount of fluid there is less. Here, again the LB results differ from MD and resemble case A.

FIG. 15.

Drop spreading radius on a cross patterned surface with Re = 24 and We = 128. (a) Case A in MD; (b) case B in MD; (c) case A in LB; (d) case B in LB. ○: wetting pattern; –: uniform wettable surface; ×: nonwetting pattern; - -: uniform nonwetting surface. Case A: ○: along axes and ×: along diagonals; Case B: ○: along diagonals and ×: along axes.

FIG. 15.

Drop spreading radius on a cross patterned surface with Re = 24 and We = 128. (a) Case A in MD; (b) case B in MD; (c) case A in LB; (d) case B in LB. ○: wetting pattern; –: uniform wettable surface; ×: nonwetting pattern; - -: uniform nonwetting surface. Case A: ○: along axes and ×: along diagonals; Case B: ○: along diagonals and ×: along axes.

Close modal

Interesting behavior occurs if we consider higher impact velocity, Re = 36, which was previously shown to produce a splash on a homogeneous nonwetting surface10 but, as discussed in Sec. III, simply deforms the drop on wettable surfaces. On the cross patterns, the result is a mixture of these behaviors: snapshots at time t = 2t* are given in Fig. 16. In case A most of the solid beneath the drop is nonwetting, to an increasing degree as the drop spreads, and the behavior in the four diagonal nonwetting directions away from the cross pattern is splash-like. Along the cross the motion is slower with a stable interface, and these segments of the drop surface act to stabilize the entire drop, somewhat mitigating the unstable splashing behavior on the diagonals. The result is an irregular pattern with holes and satellite drops, but far less disordered as that on a fully nonwetting surface. Because the drop breaks up, it does not contract at later times. In case B the surface is reversed and most of the drop spreads along a wetting surface which is stable at this Re. The result is a spreading pattern with rapidly moving fingers along the nonwetting cross axes which break up at later times and a slower and slightly irregular circular edge advancing in the diagonal directions.

FIG. 16.

Higher velocity impact in MD. Top view of the drop shape at time t = 2t* at Re = 36 and We = 228. Left: case A; right: case B. Dark region is wetting and light region is nonwetting.

FIG. 16.

Higher velocity impact in MD. Top view of the drop shape at time t = 2t* at Re = 36 and We = 228. Left: case A; right: case B. Dark region is wetting and light region is nonwetting.

Close modal

The cross-patterned surfaces provide a vivid illustration of the difference between LB and MD in fully nonwetting cases. We consider an impact at “zero velocity”: a drop is placed at rest just within interaction range of the surface (separation less than the cutoff distance 2.5σ of the Lennard-Jones interaction in MD) and allowed to move freely. The drop is drawn to the surface by the attraction from the wettable regions, but the final state depends on the pattern. Side views of the initial and final drop shapes (which stabilize after 400τ) are shown in Fig. 17. In case A the drop is attracted to the cross and nearby fluid molecules approach the solid closely. Fluid above the nonwetting solid regions would tend to float away, but the result would be a significant distortion in the drop shape and surface tension suppresses this. The resulting drop is a near-sphere, slightly flattened at the surface. In case B the cross repels nearby fluid molecules but the diagonal regions attract them, and since there is more surface along the diagonals the drop is still bound to the solid. However, as seen in the figure, a cavity opens above the cross. In contrast, LB simulations of a static drop in case B do not exhibit a cavity. Evidently, there are differences in the surface interactions between the two simulation methods: the sharp distinction between attracting and non-attracting surface atoms in MD is not entirely equivalent to a spatial variation of the boundary coefficients in the Cahn-Hilliard potential in LB.

FIG. 17.

Static drop shapes above a cross pattern. (a) Initial shape of a drop near the surface. (b) Equilibrium shape in case A. (c) Equilibrium shape in case B, viewed along the x-axis. (d) Equilibrium shape in case B, viewed along a 45° diagonal.

FIG. 17.

Static drop shapes above a cross pattern. (a) Initial shape of a drop near the surface. (b) Equilibrium shape in case A. (c) Equilibrium shape in case B, viewed along the x-axis. (d) Equilibrium shape in case B, viewed along a 45° diagonal.

Close modal

Realistic solid surfaces are not often atomically smooth, as in the examples discussed above, and a variation in height is the general case. A weak variation would presumably produce local variations in shape, while strong variations would lead to pinning effects, and a rather extensive study would be needed to explore the range of fluctuation amplitudes and correlation lengths present in different materials. Here we focus on the effects of periodic pillar arrays on spreading, since such surfaces are both practically relevant because they can exhibit super-hydrophobic behavior and theoretically interesting because of the interplay of Cassie and Wenzel states55,56 which have different degrees of liquid filling in the gaps between the pillars. In fact, some earlier MD simulations57,58 have specifically explored the transitions between the various filled states, but only for drops initially at rest.

The geometry is given in Fig. 18, and consists of square pillars of width w and height h, arranged in a square array. We vary the pillar density by changing the spacing ratio ρp = w/(w + s) while keeping w + s ≃ 0.2R0. The contact area is awkward to quantify in simulations here, since liquid may coat part of a pillar, and difficult to measure experimentally for the same reason, so we consider the projected or maximum radius Rm only.

FIG. 18.

Geometry of pillared array surfaces. (Left) Top view, where grey squares indicate the tops of the pillars. (Right) Side view, where the solid is shadowed. The pillars have width w and height h and are separated by a gap s.

FIG. 18.

Geometry of pillared array surfaces. (Left) Top view, where grey squares indicate the tops of the pillars. (Right) Side view, where the solid is shadowed. The pillars have width w and height h and are separated by a gap s.

Close modal

The pillar spacing has a strong effect on drop behavior, as shown in Fig. 19. For non-wetting pillars, when ρp = 5/8, the gap between pillars is small and the pressure in the drop is too small to overcome the capillary pressure necessary for the fluid to invade a narrow channel, giving a Cassie state of fluid advancing atop the pillars. If the pillar are more widely spaced, at ρp = 3/8, the gaps are wide enough to be invaded, resulting in a partially-filled state. Here, because the surface is nonwetting, the drop will eventually evacuate the gaps but this process delays the advance of the drop and leads to considerably smaller spreading diameters. Note that gap-filling is rather more rapid in the MD case.

FIG. 19.

Side views of spreading drops on nonwetting pillars Re = 17. First row: ρp = 5/8 in MD; second row: ρp = 5/8 in LB; third row: ρp = 3/8 in MD; bottom row: ρp = 3/8 in LB. Time from left to right: t/t* = 0.5, 1, 2, 3.

FIG. 19.

Side views of spreading drops on nonwetting pillars Re = 17. First row: ρp = 5/8 in MD; second row: ρp = 5/8 in LB; third row: ρp = 3/8 in MD; bottom row: ρp = 3/8 in LB. Time from left to right: t/t* = 0.5, 1, 2, 3.

Close modal

If instead the solid surface is wetting, the gaps are always fully invaded and produce a Wenzel state, but the time-dependent dynamics depends on the geometry. In Fig. 20, we show the drop states at the time of maximum spreading for a partially wetting surface with equilibrium contact angle 50° and two choices of pillar shape. If the pillars are shallow and widely spaced the drop first deforms to its maximum extension with a dynamic contact angle greater than 90°, filling the gaps as it spreads, and then gradually increases its contact area while decreasing the contact angle as the gap filling at the edge of the drop adjusts. However, if the pillars are deep and slightly more densely packed, the gaps are incompletely filled during spreading, and as they fill afterwards the drop withdraws to provide the necessary fluid and lowers Rm. The top view of the drop at late stage of impact shows how the surface irregularities perturb the shape of the lamella: the fluid is attracted to the pillars and the rim of the drop is flattened parallel to the axes.

FIG. 20.

Late-time evolution of drop shape on wettable pillars with θ = 50° at Re = 17 and We = 63. Top row: ρp = 2/9 and h = 0.16R0; bottom row: ρp = 3/8 and h = 0.3r0. tm is when the drop reaches its maximum deformation and the last column is the top view at t = 4.6t*.

FIG. 20.

Late-time evolution of drop shape on wettable pillars with θ = 50° at Re = 17 and We = 63. Top row: ρp = 2/9 and h = 0.16R0; bottom row: ρp = 3/8 and h = 0.3r0. tm is when the drop reaches its maximum deformation and the last column is the top view at t = 4.6t*.

Close modal

For the quantitative time dependence, first in Fig. 21, we consider the effects of varying the pillar shape, for two choices of wettability. Unsurprisingly, in both cases the spreading rates and maximum diameters on a pillared surface is less than those of a flat one. In the nonwetting case, the drop always expands and then contracts, and increasing the pillar areal density and decreasing the pillar depth both promote spreading. In the partially-wetting case the relative amount of spreading exhibits the same trend as the pillar geometry varies, but in the presence of the pillars the drops halt or contract at later times. On a flat surface this liquid would spread until the contact angle reaches the equilibrium value of 50°, at times beyond those simulated, but here the maximum spreading radius is less because liquid is “lost” in starting to fill the gaps, and the late-time decay results as the gaps fill completely.

FIG. 21.

Maximum drop spreading radius for impact on non- and partially-wetting surface for pillars of various sizes at Re = 17 and We = 63. Left: θ = 180°; right: θ = 50°. ×: ρp = 5/8; +: ρp = 3/8; △: ρp = 3/8 with shallower depth h/D0 = 0.09; −: ρp = 1 (flat surface).

FIG. 21.

Maximum drop spreading radius for impact on non- and partially-wetting surface for pillars of various sizes at Re = 17 and We = 63. Left: θ = 180°; right: θ = 50°. ×: ρp = 5/8; +: ρp = 3/8; △: ρp = 3/8 with shallower depth h/D0 = 0.09; −: ρp = 1 (flat surface).

Close modal

The interplay of pillar geometry and wettability can be illustrated by comparing the spreading curves for four contact angles in two geometries, in Fig. 22. We refer to the two pillar arrays as case C, with spacing ratio ρp = 2/9 and depth 0.16R0, and case D with ρp = 3/8 and depth 0.3R0, respectively, and the ratio of their gap volumes is about 0.6. The figure indicates a qualitative difference in behavior between the two less-wetting and the two more-wetting contact angles. On the nonwetting surfaces, drops always spread and contract, both at slower rates than in the case of a flat surface (see Fig. 11 above). The reduction in rate is due to the transition to a Wenzel state, as discussed earlier, and the fact that the reduction is similar in the two geometries even though the gap volumes differ can be attributed to incomplete filling of the gap due to the nonwetting nature of the liquid in this situation. For the two wettable surfaces, the effect of the pillars is again to halt the advance after the initial spreading, at times around 2t*, but the subsequent behavior in the two geometries differs. Here, the liquid will fill the gaps completely, but more is required to fill the higher volume in case D, so that in case C there is additional liquid available for the drop to continue to spread to reach its final contact angle.

FIG. 22.

Drop spreading radius on pillars of various wettabilities with Re = 17 and We = 63. (a) ρp = 2/9 and depth 0.16R0; (b) ρp = 3/8 and depth 0.3R0. +: θ = 180°; ○: θ = 96°; △: θ = 50°; □: θ = 0°.

FIG. 22.

Drop spreading radius on pillars of various wettabilities with Re = 17 and We = 63. (a) ρp = 2/9 and depth 0.16R0; (b) ρp = 3/8 and depth 0.3R0. +: θ = 180°; ○: θ = 96°; △: θ = 50°; □: θ = 0°.

Close modal

A different form of deviation from a flat surface is curvature, which may be though of as a model of large scale roughness appropriate to a drop which is smaller than the characteristic size of surface features. There are again numerous distinct impact configurations available, depending on the local value of the curvature and the orientation of the impact velocity with respect to the surface, but we focus on the simplest cases of normal impact on a concave or convex but otherwise homogeneous surface. Examples of both cases are shown in Fig. 23 for two choices of wettability and a radius of curvature of 4R0.

FIG. 23.

Side view of spreading drops on walls with radius of curvature 4R0. The solid and dotted lines lines are the drop and wall shapes, respectively, with concave in dark gray (red) and convex in light gray (blue). Top row: θ = 180°; bottom row: θ = 50°. Time from left to right: t/t* = 1, 2, 5.

FIG. 23.

Side view of spreading drops on walls with radius of curvature 4R0. The solid and dotted lines lines are the drop and wall shapes, respectively, with concave in dark gray (red) and convex in light gray (blue). Top row: θ = 180°; bottom row: θ = 50°. Time from left to right: t/t* = 1, 2, 5.

Close modal

Generally the spreading drop attempts to move along the surface and its bottom follows the local curvature for both wettabilities, even after the drop floats off the nonwettable concave surface. Presumably this is an an effect of inertia: the drop molecules' initial downward momentum is first redirected laterally by the impact, and a concave shape acts to redirect this lateral momentum upwards and towards the drop axis, carrying the drop upwards and slightly inwards. In the convex case, instead, some of the downward component of the initial momentum persists and pushed the drop downwards along the surface. The drop then spreads further along the surface and the lamella thins in comparison to the concave case. The quantitative results for the thickness are given in Fig. 24, which indicates that convex and nonwettable surfaces have thinner lamella than concave wettable ones.

FIG. 24.

Drop thickness at the center after impact on curved surfaces at We = 128 and Re = 24. Left: θ = 180°; right: θ = 50°. Open (blue) symbols refer to a convex surface and filled (red) markers are for concave and the solid lines refer to a flat surface. Radius of curvature: ◊ or ♦ ✦: R = ±2D0; ○ or •: R = ±6D0; △ or ▲: R = ±10D0; −: R = ∞.

FIG. 24.

Drop thickness at the center after impact on curved surfaces at We = 128 and Re = 24. Left: θ = 180°; right: θ = 50°. Open (blue) symbols refer to a convex surface and filled (red) markers are for concave and the solid lines refer to a flat surface. Radius of curvature: ◊ or ♦ ✦: R = ±2D0; ○ or •: R = ±6D0; △ or ▲: R = ±10D0; −: R = ∞.

Close modal

The influence of the value of the curvature on spreading is shown in Fig. 25, for three different wettabilities (the completely wetting case is similar to that at 50°). The spreading radius in this situation is taken to be the circumferential distance along the surface. If we take the sign of the curvature to be negative for concave and positive for convex, the overall qualitative effect is that both the maximum spreading radius and the time required to reach it increase with curvature. More precisely, when these quantities are plotted directly in Fig. 26, there is a very strong variation with curvature on nonwetting surfaces, for the reasons discussed above. In the wetting case the liquid tends to wrap itself around the solid whatever its shape and spread until the impact kinetic energy and momentum are dissipated, almost independently of the curvature.

FIG. 25.

Drop spreading on surfaces with different curvature at Re = 24 and We = 128. (a) θ = 180°; (b) θ = 96°; (c)θ = 50°. Open (blue) markers are convex surfaces, filled (red) markers are concave surfaces and the solid line is for a flat surface. The radii of curvature shown are ◊ or ✦ ♦: R = ±4R0: ○ or •: R = ±12R0; △ or ▲: R = ±20R0.

FIG. 25.

Drop spreading on surfaces with different curvature at Re = 24 and We = 128. (a) θ = 180°; (b) θ = 96°; (c)θ = 50°. Open (blue) markers are convex surfaces, filled (red) markers are concave surfaces and the solid line is for a flat surface. The radii of curvature shown are ◊ or ✦ ♦: R = ±4R0: ○ or •: R = ±12R0; △ or ▲: R = ±20R0.

Close modal
FIG. 26.

(Left) Maximum spreading factor ξmax on surface of curvature R0/R; (Right) Time to reach maximum deformation tm/t* on surface of curvature R0/R. Wettability: △: θ = 180°; □: θ = 180°; ◊: θ = 50°. tm for θ = 50° is not shown.

FIG. 26.

(Left) Maximum spreading factor ξmax on surface of curvature R0/R; (Right) Time to reach maximum deformation tm/t* on surface of curvature R0/R. Wettability: △: θ = 180°; □: θ = 180°; ◊: θ = 50°. tm for θ = 50° is not shown.

Close modal

In this paper we have used molecular dynamics and lattice Boltzmann simulation techniques to simulate drop impact on surfaces with various shapes and textures for a range of moderate impact velocities ( 10 ⩽ Re ⩽ 36 and 20 ⩽ We ⩽ 300) and fixed Ohnesorge number Oh = 0.4774. Lower velocities lead to smaller initial deformations and dynamics dominated by quasi-static spreading considerations. Higher velocities lead to drop splashing, for which neither method is appropriate. Aside from illustrating the evolving drop shapes, we quantify the behavior in terms of the time dependence of the drop height and spreading radius. Since MD operates at atomic length and time scales while LB is constructed to solve the Navier-Stokes equations at continuum scales, a comparison of the results of the two methods is of some interest. We find reasonable agreement between the methods for slow impact on wettable surfaces, and some deviation otherwise. There are two sources of disagreement: compressibility in the MD case and boundary effects in LB. In order to increase the Reynolds number to higher values for a nanometer-sized drop the impact velocity exceeds the thermal velocity of the fluid and begins to approach the sound speed. The observed variation in volume at the highest impact velocities discussed in this paper is around 5% at Ma ≈ 0.3. In principle the LB calculations can be modified to finite Ma, but convergence issues arise at the relevant high velocities. The second distinction between the methods is in the treatment of solid boundaries. Both MD and LB work well for surfaces which are at least partially wettable, where a no-slip boundary condition is appropriate at low shear rates. Otherwise, MD predicts systematically increasing slip lengths and, for very low wettability, only tenuous contact between the drop and the surface. In LB the surface interaction is implemented in an entirely different way, and the slip velocity is independent of wettability during the spreading stage.

Aside from a comparison of two methods, we have explored the behavior of drops impacting surfaces which are not atomically smooth and homogeneous. To study the effects of the variability of surface wetting properties, we considered an example studied previously involving a surface with a cross pattern of one wettability embedded in a differently wettable background. The fact that spreading is more rapid as wettability decreases leads to final drop shapes with various indentations and protuberances, and different degrees of spreading. In the case of surface structural heterogeneity, we focused on periodic arrays of pillars, a configuration relevant to super-hydrophobic behavior. Here the filling of the gaps between pillars leads to lower areal coverage by the drop, related to transitions between Cassie and Wenzel states. Finally, we considered surfaces with various curvatures which reflect larger scale rough features. Curvature has a pronounced effect in the nonwetting case, and leads to spreading distances and times which increase with curvature, a thinner spreading lamella and more of a tendency for an impacting drop to bounce, although the effect is small in the wetting case.

The results in this paper indicate a partial overlap in the results of the two computational methods, along with distinct limitations on the range of operating conditions where this agreement holds. The MD method does not have many internal “computational” parameters to adjust: once the potential is specified, quantities such as the time step are fixed by numerical accuracy and stability requirements, the material properties and transport coefficients are determined, and the numbers of the various atomic species present are in direct correspondence to the physical size of the system being studied. The LB method is, like the Navier-Stokes equations themselves, not restricted to a particular system scale, but it does involve a number of internal variables without a direct physical counterpart, such as the Cahn-Hilliard potential and the mobility, which allow the material and dynamic properties to be adjusted to some degree. In this paper, the particular LB algorithm has allowed some improvements along these lines, which might be regarded as calibrating the method using the results of (external, numerical) experiments, but a robust computational tool spanning the full range of hydrodynamic phenomena needs further development.

We do not have a general answer to the question of which method is “correct.” The MD calculations presumably capture the properties of nanoscale drops of a generic viscous liquid in high velocity impacts on the modeled solid substrate, but do not necessarily apply when extrapolated to the regime where experiments are carried out. The fact that the LB results are somewhat different means that the implementation of the method used here may not be entirely accurate in the nano domain. In fact this method is explicitly meant to reproduce the (macroscopic) Navier-Stokes equation and indeed it is generally in better agreement with the corresponding experiments.

We thank Kevin Connington for discussions. This work is supported in part by a NSF PREM Grant No. DMR0934206.

1.
H.
Dong
,
W. W.
Carr
,
D. G.
Bucknall
, and
J. F.
Morris
, “
Temporally-resolved inkjet drop impaction on surfaces
,”
AIChE J.
53
,
2606
2617
(
2007
).
2.
T.
Bennett
and
D.
Poulikakos
, “
Splat-quench solidification: estimating the maximum spreading of a droplet impacting a solid surface
,”
J. Mater. Sci.
28
,
963
970
(
1993
).
3.
L.
Xu
,
W. W.
Zhang
, and
S. R.
Nagel
, “
Drop splashing on a dry smooth surface
,”
Phys. Rev. Lett.
94
,
184505
(
2005
).
4.
A.
Theodorakakos
,
T.
Ous
,
M.
Gavaises
,
J. M.
Nouri
,
N.
Nikolopoulos
, and
H.
Yanagihara
, “
Dynamics of water droplets detached from porous surfaces of relevance to PEM fuel cells
,”
J. Colloid Interface Sci.
300
,
673
687
(
2006
).
5.
C.
Antonini
,
M.
Innocenti
,
T.
Horn
,
M.
Marengo
, and
A.
Amirfazli
, “
Understanding the effect of superhydrophobic coatings on energy reduction in anti-icing systems
,”
Cold Reg. Sci. Technol.
67
,
58
67
(
2011
).
6.
P.
Tsai
,
R.
van der Veen
,
M.
van de Raa
, and
D.
Lohse
, “
How micropatterns and air pressure affect splashing on surfaces
,”
Langmuir
26
,
16090
16095
(
2010
).
7.
J. B.
Lee
and
S. H.
Lee
, “
Dynamic wetting and spreading characteristics of a liquid droplet impinging on hydrophobic textured surfaces
,”
Langmuir
27
,
6565
6573
(
2011
).
8.
A. L.
Yarin
, “
Drop impact dynamics: Splashing, spreading, receding, bouncing…
,”
Annu. Rev. Fluid Mech.
38
,
159
192
(
2006
).
9.
R.
Rioboo
,
C.
Tropea
, and
M.
Marengo
, “
Outcomes from a drop impact on solid surfaces
,”
Atomization Spray
11
,
2
(
2001
).
10.
J.
Koplik
and
R.
Zhang
, “
Nanodrop impact solid surfaces
,”
Phys. Fluids
25
,
022003
(
2013
).
11.
W.
Chen
,
R.
Zhang
, and
J.
Koplik
, “
Velocity slip on curved surface
,”
Phys. Rev. E
89
,
023005
(
2014
).
12.
H.
Marmanis
and
S. T.
Thoroddsen
, “
Scaling of the fingering pattern of an impacting drop
,”
Phys. Fluids
8
,
1344
1346
(
1996
).
13.
R.
Dhiman
and
S.
Chandra
, “
Rupture of thin films formed during droplet impact
,”
Proc. R. Soc. London, Ser. A
466
,
1229
1245
(
2010
).
14.
L.
Courbin
,
J. C.
Bird
,
A.
Belmonte
, and
H. A.
Stone
, “
‘Black hole' nucleation in a splash of milk
,”
Phys. Fluids
20
,
091106
(
2008
).
15.
J. G.
Leidenfrost
,
De Aquae Communis Nonnul lis Qualitatibus Tractatus
(
Duisburg
,
1756
)
[English translation by
C.
Wares
,
Int. J. Heat Mass Transfer
9
,
1153
(
1966
)].
16.
D.
Quéré
, “
Leidenfrost dynamics
,”
Annu. Rev. Fluid Mech.
45
,
197
215
(
2013
).
17.
F.
Celestini
,
T.
Frisch
, and
Y.
Pomeau
, “
Take off of small Leidenfrost droplets
,”
Phys. Rev. Lett.
109
,
034501
(
2012
).
18.
J.
Horbach
and
S.
Succi
, “
Lattice Boltzmann versus molecular dynamics simulation of nanoscale hydrodynamic flows
,”
Phys. Rev. Lett.
96
,
224503
(
2006
).
19.
L. W.
Schwartz
and
R. R.
Eley
, “
Simulation of droplet motion on low-energy and heterogeneous surfaces
,”
J. Colloid Interface Sci.
202
,
173
188
(
1998
).
20.
L. W.
Schwartz
, “
Hysteretic effects in droplet motions on heterogeneous substrates: Direct numerical simulation
,”
Langmuir
14
,
3440
3453
(
1998
).
21.
V.
Vaikuntanathan
,
R.
Kannan
, and
D.
Sivakumar
, “
Impact of water drops onto the junction of a hydrophobic texture and a hydrophilic smooth surface
,”
Colloids Surf., A
369
,
65
74
(
2010
).
22.
L.
Cheng
, “
Dynamic spreading of drops impacting onto a solid surface
,”
Ind. Eng. Chem. Process Des. Dev.
16
,
192
(
1977
).
23.
U.
Mock
,
T.
Michel
,
C.
Tropea
,
I.
Roisman
, and
J.
Rühe
, “
Drop impact on chemically structured arrays
,”
J. Phys. Condens. Matter
17
,
S595
S605
(
2005
).
24.
M.
Lee
,
Y. S.
Chang
, and
H.-Y.
Kim
, “
Drop impact on microwetting patterned surfaces
,”
Phys. Fluids
22
,
072101
(
2010
).
25.
S.
Kim
,
M.-W.
Moon
, and
H.-Y.
Kim
, “
Drop impact on super-wettability-contrast annular patterns
,”
J. Fluid Mech.
730
,
328
342
(
2013
).
26.
H.
Katsuragi
, “
Morphology scaling of drop impact onto a granular layer
,”
Phys. Rev. Lett.
104
,
218001
(
2010
).
27.
P.
Tsai
,
M. W.
Hendrix
,
R. M.
Dijkstra
,
L.
Shui
, and
D.
Lohse
, “
Microscopic structure influencing macroscopic splash at high Weber number
,”
Soft Matter
7
,
11325
(
2011
).
28.
M.
Sbragaglia
,
A. M.
Peters
,
C.
Pirat
,
B. M.
Borkent
,
R. G. H.
Lammertink
,
M.
Wessling
, and
D.
Lohse
, “
Spontaneous breakdown of superhydrophobicity
,”
Phys. Rev. Lett.
99
,
156001
(
2007
).
29.
G. O.
Berim
and
E.
Ruckenstein
, “
Nanodrop on a nanorough solid surface: Density functional theory considerations
,”
J. Chem. Phys.
129
,
014708
(
2008
).
30.
G. O.
Berim
and
E.
Ruckenstein
, “
Nanodrop on a nanorough hydrophilic solid surface: Contact angle dependence on the size, arrangement, and composition of the pillars
,”
J. Colloid Interface Sci.
359
,
304
310
(
2011
).
31.
J.-L.
Barrat
and
L.
Bocquest
, “
Large slip effect at a nonwetting fluid-solid interface
,”
Phys. Rev. Lett.
82
,
4671
(
1999
).
32.
G. S.
Grest
and
K.
Kremer
, “
Molecular dynamics simulation for polymers in the presence of a heat bath
,”
Phys. Rev. A
33
,
3628
(R) (
1986
).
33.
W.
Bouwhuis
,
R. C. A.
van der Veen
,
T.
Tran
,
D. I.
Keij
,
K. G.
Winkels
,
I. R.
Peters
,
D.
van der Meer
,
C.
Sun
,
J. H.
Snoeijer
, and
D.
Lohse
, “
Maximal air bubble entrainment at liquid-drop impact
,”
Phys. Rev. Lett.
109
,
264501
1
264501
4
(
2012
).
34.
A.
Dupuis
and
J. M.
Yeomans
, “
Modeling droplets on superhydrophobic surfaces: equilibrium states and transitions
,”
Langmuir
21
,
2624
(
2005
).
35.
H.
Kusumaatmaja
,
M. L.
Blow
,
A.
Dupuis
, and
J. M.
Yeomans
, “
The collapse transition on superhydrophobic surfaces
,”
Europhys. Lett.
81
,
36003
(
2008
).
36.
R. J.
Vrancken
,
H.
Kusumaatmaja
,
K.
Hermans
,
A. M.
Prenen
,
O.
Pierre-Louis
,
C. W. M.
Bastiaansen
, and
D. J.
Broer
, “
Fully reversible transition from Wenzel to Cassie-Baxter states on corrugated superhydrophobic surfaces
,”
Langmuir
26
,
3335
(
2010
).
37.
J.
Zhang
and
D. Y.
Kwok
, “
Contact line and contact angle dynamics in superhydrophobic channels
,”
Langmuir
22
,
4998
(
2006
).
38.
J. J.
Huang
,
C.
Shu
, and
Y. T.
Chew
, “
Lattice Boltzmann study of droplet motion inside a grooved channel
,”
Phys. Fluids
21
,
022103
(
2009
).
39.
N.
Moradi
,
F.
Varnik
, and
I.
Steinbach
, “
Roughness-gradient-induced spontaneous motion of droplets on hydrophobic surfaces: a lattice Boltzmann study
,”
Europhys. Lett.
89
,
26006
(
2010
).
40.
K.
Connington
and
T.
Lee
, “
Lattice Boltzmann simulations of forced wetting transitions of drops on superhydrophobic surfaces
,”
J. Comput. Phys.
250
,
601
615
(
2013
).
41.
T.
Lee
and
L.
Liu
, “
Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces
,”
J. Comput. Phys.
229
,
8045
8063
(
2010
).
42.
T.
Lee
, “
Effects of incompressibility on the elimination of parasitic currents in the lattice Boltzmann equation method for binary fluids
,”
Comput. Math. Appl.
58
,
987
994
(
2009
).
43.
L.
Barbieri
,
E.
Wagner
, and
P.
Hoffmann
, “
Water wetting transition parameters of perfluorinated substances with periodically distributed flat-top microscale obstacles
,”
Langmuir
23
,
1723
(
2007
).
44.
L.
Liu
and
T.
Lee
, “
Wall free energy based polynomial boundary conditions for non-ideal gas lattice boltzmann equation
,”
Int. J. Mod. Phys. C
20
(
11
),
1749
1768
(
2009
).
45.
J. W.
Cahn
, “
Critical point wetting
,”
J. Chem. Phys.
66
,
3667
(
1977
).
46.
K. G.
Winkels
,
J. H.
Weijs
,
A.
Eddi
, and
J. H.
Snoeijer
, “
Initial spreading of low-viscosity drops on partially wetting surfaces
,”
Phys. Rev. E
85
,
055301
1
055301
4
(
2012
).
47.
C.
Clanet
,
C.
Béguin
,
D.
Richard
, and
D.
Quéré
, “
Maximal deformation of an impacting drop
,”
J. Fluid Mech.
517
,
199
(
2004
).
48.
P.
Yue
,
C.
Zhou
, and
J. J.
Feng
, “
Sharp-interface limit of the Cahn-Hilliard model for moving contact lines
,”
J. Fluid Mech.
645
,
279
(
2010
).
49.
R. G.
Cox
, “
The dynamics of the spreading of liquids on a soild surface. Part I. Viscous flow
,”
J. Fluid Mech.
168
,
169
194
(
1986
).
50.
B. L.
Scheller
and
D. W.
Bousfield
, “
Newtonian drop impact with a solid surface
,”
AICHE J.
41
,
1357
1367
(
1995
).
51.
I. V.
Roisman
, “
Inertia dominated drop collisions. II. An analytical solution of the Navier-Stokes equations for a spreading viscous film
,”
Phys. Fluids
21
,
052104
(
2009
).
52.
M.
Pasandideh-Fard
,
Y.
Qiao
,
S.
Chandra
, and
J.
Mostaghimi
, “
Capillary effects during droplet impact on a solid surface
,”
Phys. Fluids
8
,
650
(
1996
).
53.
C. W.
Visser
,
Y.
Tagawa
,
C.
Sun
, and
D.
Lohse
, “
Microdroplet impact at very high velocity
,”
Soft Matter
8
,
10732
10737
(
2012
).
54.
S.
Chandra
and
C. T.
Avedisian
, “
On the collision of a droplet with a solid surface
,”
Proc. R. Soc. London, Ser. A
432
,
13
41
(
1991
).
55.
R. N.
Wenzel
, “
Resistance of solid surfaces to wetting by water
,”
Ind. Eng. Chem.
28
,
988
994
(
1936
).
56.
A. B. D.
Cassie
and
S.
Baxter
, “
Wettability of porous surfaces
,”
Trans. Faraday Soc.
40
,
546
551
(
1944
).
57.
T.
Koishi
,
K.
Yasuoka
,
S.
Fujikawa
,
T.
Ebisuzaki
, and
X. C.
Zeng
, “
Coexistence and transition between Cassie and Wenzel state on pillared hydrophobic surface
,”
Proc. Natl. Acad. Sci. U.S.A.
106
,
8435
8440
(
2009
).
58.
A.
Giacomello
,
S.
Meloni
,
M.
Chinappi
, and
C. M.
Casciola
, “
Cassie-Baxter and Wenzel states on a nanostructured surface: Phase diagram, metastabilities, and transition mechanism by atomistic free energy calculations
,”
Langmuir
28
,
10764
10772
(
2012
).