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.

## I. INTRODUCTION

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 effect^{15,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 liquid^{18} 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. Cheng^{22} 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 Katsuragi^{26} 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.

## II. SIMULATION METHODS

### A. Molecular dynamics

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,

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 − 2*c*ρ_{s}/ρ_{l},^{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.

c | 0 | 0.4 | 0.67 | 0.8 | 0.9 | 1.0 |

θ | 180° | 134° | 96° | 73° | 50° | 0° |

c | 0 | 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ε/*k*_{B}, 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 *v*_{0}, 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 *R*_{m}, and the radius of the contact area *R*_{c}. 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 *R*_{0}.

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 *R*_{m} are now well defined, there is some further ambiguity in the measurement of *R*_{c} 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 *b*_{c} 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*) < *b*_{c}.

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* = 2*R*_{0}. 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.

### B. Lattice Boltzmann method

The LBM was first used by Dupuis and Yeomans^{34} 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 Lee^{40} presented an extension of the method of Lee and Liu^{41,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 Liu^{41,42} using a scalar order parameter *C*, which satisfies the convective Cahn-Hilliard equation ∂_{t}*C* + ∇ · (**u***C*) = *M*∇^{2}μ, 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

involving a bulk free energy *E*_{0}(*C*) = β*C*^{2}(*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 *C*_{s}. 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

_{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 *D*3*Q*27 lattice. The boundaries are all symmetric except at the solid surface, where the wall boundary condition^{44} is imposed. The interface thickness and initial drop radius are ζ = 5 and *R*_{0} = 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:

where *u*_{0} is the drop impact speed, *D*_{0} 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 *D*2*Q*9 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*/*D*_{0} 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.

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, *D*_{max} ∼ *D*_{0}*We*^{1/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.

## III. HOMOGENEOUS SURFACES

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** = *R*_{0}/*v*_{0}, 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.

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 *R*_{0} and time is scaled (inertially) by the free-space transit time across a radius, *t** = *R*_{0}/*v*_{0}. 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.

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 *R*_{m}, corresponding to the drop size as projected from above, and the contact radius of liquid on the surface *R*_{c}, 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 *R*_{0} 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.

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 surfaces^{10} 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 shown^{48} 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 Cox^{49}). 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.

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.

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.

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 *R*_{c} varies smoothly and time-dependence of *R*_{m} 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-3*t** 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 *t*^{1/2} growth law for the contact radius. Our results are consistent with this behavior.

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 ξ_{max} ∼ *We*^{1/4},^{47} where ξ_{max} is the maximum spreading diameter during impact scaled by the initial drop diameter *D*_{0}. An empirical extrapolation from experiments by Scheller and Bousfield,^{50} a combination of data and theory by Roisman^{51} 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 *We*^{0.49} for *D*_{m} and *We*^{0.56} for *D*_{c} (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.

## IV. CROSS PATTERNS

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 *R*_{0}/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.

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* ⩽ 2*t** 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.

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 3*t**, 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.

Interesting behavior occurs if we consider higher impact velocity, *Re* = 36, which was previously shown to produce a splash on a homogeneous nonwetting surface^{10} 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* = 2*t** 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.

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.

## V. PERIODIC PILLAR ARRAYS

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 states^{55,56} which have different degrees of liquid filling in the gaps between the pillars. In fact, some earlier MD simulations^{57,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.2*R*_{0}. 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 *R*_{m} only.

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.

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 *R*_{m}. 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.

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.

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.16*R*_{0}, and case D with ρ_{p} = 3/8 and depth 0.3*R*_{0}, 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 2*t**, 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.

## VI. CURVED SURFACES

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 4*R*_{0}.

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.

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.

## VII. CONCLUSION

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.

## ACKNOWLEDGMENTS

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