Numerical Simulations of the Impact and Spreading of a Particulate Drop on a Solid Substrate

Abstract

We present two-dimensional numerical simulations of the impact and spreading of a droplet containing a number of small particles on a flat solid surface, just after hitting the solid surface, to understand particle effects on spreading dynamics of a particle-laden droplet for the application to the industrial inkjet printing process. The Navier-Stokes equation is solved by a finite-element-based computational scheme that employs the level-set method for the accurate interface description between the drop fluid and air and a fictitious domain method for suspended particles to account for full hydrodynamic interaction. Focusing on the particle effect on droplet spreading and recoil behaviors, we report that suspended particles suppress the droplet oscillation and deformation, by investigating the drop deformations for various Reynolds numbers. This suppressed oscillatory behavior of the particulate droplet has been interpreted with the enhanced energy dissipation due to the presence of particles.

1. Introduction

In this work, we present numerical simulations for the impact and spreading of a particulate droplet, containing a large number of solid particles, on a flat substrate and focus on particle effects on droplet dynamics such as spreading and recoil behaviors. Understanding such a behavior are particularly important in the industrial inkjet printing technology, which provides an alternative to the conventional lithography process for the purpose of accurate particle delivery to a solid surface. The present work is an extension of the authors’ previous study on a particulate droplet deformation/spread under the capillary force [1] to the inertial spreading problem by solving the Navier-Stokes equation during the droplet impact.

Numerous studies on droplet impact and spreading on a solid surface have been reported. Fukai and coworkers [2, 3] reported experimental and theoretical studies on the deformation of a droplet—the splat radius and height—impinging on a flat surface at a high Reynolds number (Re≈4500) in an axisymmetric moving mesh. Šikalo et al. [4] performed experiments on drop impact and spreading and investigated the effective contact line and similar experiments by Park et al. [5] using a spherical cap improved predictions for flows with low Reynolds number, low Weber number, and a large contact angle. Pasandideh-Fard et al. [6] and Gunjal et al. [7] performed numerical simulations on the droplet impact using a volume-of-fluid (VOF) method and Khatavkar [8] and Khatavkar et al. [9] employed a diffuse-interface method to study the capillary and low inertia spreading of a micro droplet on a flat and patterned surface at low Reynolds numbers. Recently, Lee and Son [10] presented droplet impact and coalescence to understand microline patterning process with the level-set method. In most previous studies, however, both the experiments and simulations focused on drop deformation with a homogeneous fluid, whose application to the particulate fluid drop is limited due to the presence of the hydrodynamic interaction between particles and the fluid. Although interest in the particle effects as they pertain to drop deformation has increased, only a few experimental studies have explored the particulate droplet behavior. Furbank and Morris [11] observed the drop formation characteristics and Nicholas [12] studied the effect of particles on the spreading behavior, finding that a particulate drop spread less than a pure drop. Recently, the authors performed 2D direct simulations of low speed capillary spreading of a particulate droplet considering hydrodynamic interaction of particles and equilibrium contact angle [1].

In this work, we develop a direct numerical simulation technique for particulate drop spreading with the impact on a solid substrate. This work is restricted to two-dimensional simulations, owing to the complexity in the droplet impact and spread phenomena of the particulate drop. We present numerical results for the evolution of the drop shape, the particle distribution and movement, and the effects of the particles on the dynamics of droplet spreading. The term “direct” is used here because hydrodynamic and interparticle interactions are fully considered at the individual particle level. For this purpose, we employ the direct simulation scheme of Hwang et al. [13] similar to the distributed Lagrangian-multiplier (DLM) method [14]. We use the second-order Adams-Bashforth/Crank-Nicholson method (AB2/CN) to solve the Navier-Stokes equation and the level-set method to track the drop-air interface, dealing with a free-surface problem in a fixed mesh. The interfacial tension is assigned via the surface stress, known as the continuous surface stress (CSS) formulation and a virtual stress tensor is introduced at the contact point to exert force toward the equilibrium contact angle. Using the proposed numerical method, we study the impact and spreading behaviors of a particulate drop on a flat surface and discuss the effect of such particles by investigating particle movements, changes in the drop spreading and oscillation behaviors, and the kinetic energy.

2. Modeling

In the present work, we consider the impact and spreading of a 2D particulate liquid drop on a solid substrate and, as described in Figure 1, an initially circular drop with radius 𝑅^𝑖_𝑑 is in contact with a solid surface. An instantaneous artificial gravity is applied for a short time period to generate the initial velocity distribution within a drop as well as in the air domain, satisfying the incompressibility condition of the entire domain. The resultant velocity at the drop center, generated in this way, is regarded as the characteristic drop impact velocity 𝑉₀. Therefore the drop containing a number of randomly distributed rigid circular particles undergoes impact, spreading, and recoil behaviors through the inertial force, the elastic interfacial force, and the capillary force and subsequently reaches an equilibrium state. The spreading of a droplet in 2D is a free-surface problem with a moving contact point that requires a proper representation of various interfaces, such as the air-liquid and liquid-solid interfaces, to account for the interfacial tension and wetting effects [15]. In addition, particles inside the drop experience complex flow fields during this process, and particle motions are determined by hydrodynamic interaction between particles and the surrounding drop fluid.

Figure 1 

A schematic representation of a particulate droplet spreading. The initial drop contacts the solid surface and the particles are distributed inside the drop. The drop spreads with inertial force and then spreads to an equilibrium state.

The problem is described in Figure 2 along with the coordinate system. The computational domain of the size 𝐿×𝐻 is denoted by Ω. The four boundaries of the domain are denoted by Γ_𝑖(𝑖=1,2,3,4). The region occupied by a droplet is denoted by 𝐷(𝑡) with the interface Φ between the drop and the air. The symbol 𝐧 is the outward unit normal vector on the drop interface. Particles are denoted by 𝑃_𝑖(𝑡) (i=1,…,N) and 𝑁 is the number of particles. We use a symbol 𝑃(𝑡) for ⋃^𝑁_𝑖=1𝑃_𝑖(𝑡), a collective region occupied by the particles at a certain time 𝑡. For a particle 𝑃_𝑖(𝑡), 𝐗_𝑖=(𝑋_𝑖,𝑌_𝑖),𝐔_𝑖=(𝑈_𝑖,𝑉_𝑖), 𝝎_𝑖=𝜔_𝑖𝐤,Θ_𝑖=Θ_𝑖𝐤, and 𝑅_𝑖 are used for the locations of the particle center, the translational velocity, the angular velocity, the angular rotation, and the radius, respectively. The vector 𝐤 is the unit vector normal to the plane.

Figure 2 

The computational domain with the droplet interface and particles inside in this study.

The set of governing equations for the region 𝐷(𝑡) occupied by a drop at a certain time 𝑡 is given by 𝜌_𝑑𝑑𝐮𝑑𝑡=𝜌_𝑑(𝑡),(𝑡),𝐠+∇⋅𝝈in𝐷∇⋅𝐮=0,in𝐷𝝈=−𝑝𝐈+2𝜂_𝑑𝐃,in𝐷(𝑡). (1) Equations are equations for the momentum balance, the continuity, and the constitutive relationship of a Newtonian fluid and the symbols 𝜌_𝑑, 𝐮, 𝐠, 𝝈, 𝑝, 𝐈, 𝜂_𝑑, and 𝐃 are the drop density, the velocity, the gravity, the stress, the pressure, the identity tensor, the drop viscosity, and the rate of deformation tensor, respectively. Accepting the Eulerian approach, we replace the air with a fictitious fluid having a viscosity of one hundredth of the droplet viscosity [16]. Then we have exactly the same set of equations in the air domain (the fictitious fluid) for the momentum balance, the continuity, and the constitutive equations as those of the droplet fluid and the air region is identified by the air viscosity 𝜂_𝑎 and the air density 𝜌_𝑎. Note that the subscript “𝑎” denotes the air domain. In this sense and for convenience in future use, one can introduce quantities 𝜌^∗ and 𝜂^∗ such that 𝜌^∗=𝜌_𝑑 and 𝜂^∗=𝜂_𝑑 in the drop domain and 𝜌^∗=𝜌_𝑎 and 𝜂^∗=𝜂_𝑎 in the air domain. This Eulerian approach facilitates the use of single mesh for the entire computations with the presence of particles and the fluid interface.

As mentioned earlier, to apply the inertial force to a drop, we employ an instantaneous artificial gravity, and it is applied to an initially stationary drop only for a short time period to accelerate the drop and to achieve a specific initial velocity 𝑉₀. For this purpose, we replace the gravity 𝐠 with the instantaneous artificial gravity 𝐠 in the momentum equation as follows: 𝜌^∗𝑑𝐮𝑑𝑡=𝜌^∗𝐠+∇⋅𝝈. (2) As for the interface condition between the drop and the air, the following stress balance should be satisfied in the absence of interfacial tension gradient: 𝐧⋅𝝈_𝑎−𝝈_𝑑=𝜍𝜅𝐧,onΦ. (3) The variables 𝜍 and 𝜅 denote the interfacial tension and the local mean curvature, respectively, on the interface and 𝐧 is the local unit outward normal vector to the interface. To deal with the interface condition (3) in the Eulerian frame, we employ the level-set method and in this description the interface is not explicitly tracked, but is conveniently described by the level-set function 𝜑. The zero-level set indicates the interface and it takes a positive value within the air region and a negative value inside the drop. In fact, the level-set function is a signed distance function such that the magnitude 𝜑(𝐱) is determined by the shortest distance from a point 𝐱 to the interface: that is, ∇𝜑=1,𝜑∣_Φ=0. (4) This level-set function is advected passively under a given flow field and therefore its evolution in the Eulerian frame can be expressed as 𝜕𝜑𝜕𝑡+𝐮⋅∇𝜑=0,𝜑∣_𝑡=0=𝜑₀,inΩ. (5) If the level-set function is seriously distorted as time evolves, the level-set function must be redistributed to satisfy the distance function property (4) and we employ a direct geometric redistancing scheme developed by Kim and Hwang [17] for this purpose. This level-set description allows a simple treatment of the interfacial condition (3) by introducing the smeared-out interfacial stress 𝝉_𝑠 with the delta function 𝛿(𝜑) based on the level-set function, which is called the continuous surface stress (CSS) formulation [18]: 𝝉_𝑠=𝜍(𝐈−𝐧𝐧)𝛿𝜑, (6) and it can be added to the constitutive equation as an additional stress 𝐓=𝝈+𝝉_𝑠=−𝑝𝐈+2𝜂^∗𝐃+𝝉_𝑠. (7) where 𝐓 is the total stress and (7) defines the constitutive relation over the entire domain, including the droplet/air interface and, as we will see later, the interior of suspended particles. In addition due to our unified descriptions for drop and air domain, the momentum balance for the entire domain can be written simply as 𝜌^∗𝑑𝐮𝑑𝑡=𝜌^∗𝐠+∇⋅𝐓inΩ. (8)

In this problem, there are three driving forces in the droplet dynamics: the impact of droplet, elasticity associated with the interfacial tension, and the capillary force. Among them, the capillary force is the microscale action on the contact point with a solid wall, where complicated microscopic phenomena are present, leading to in an equilibrium contact angle 𝜃_eq macroscopically. In the present work, only the static equilibrium is considered following the work of Bellet [19], where an additional force at the contact point is introduced, if the current contact angle is different from the static equilibrium contact angle 𝜃_eq, to generate the motion toward static equilibrium. To combine it with the CSS formulation, we artificially enforce the interfacial stress tensor 𝝉_𝑠 near the contact point (within around an element distance) such that 𝝉_𝑠=𝜍𝐈−𝐧_eq𝐧_eq𝛿𝜑,𝐧_eq=cos𝜃_eq𝐧_𝑤+sin𝜃_eq𝐭_𝑤, (9) using the equilibrium normal vector 𝐧_eq, the normal 𝐧_𝑤, and tangential 𝐭_𝑤, vectors to the solid wall (Figure 3).

Figure 3 

Modification of the normal vector at the contact point.

Now we consider the description of freely suspended particles. Following the work by Hwang et al. [13], a circular particle is modeled as a rigid ring filled with the same fluid used in the drop fluid domain. The rigid body condition is imposed on the particle boundary only and it provides a plausible description for particles in fluid in low Reynolds number limit. This fictitious domain approach for the particle description is particularly useful in this work, given that it facilitates the use of a single Eulerian mesh for unified treatments of both the free-surface with interfacial tension and the solid/liquid interaction. The rigid-body condition on the particle boundary states that the velocity of a fluid particle attached onto the particle boundary is identical to that of the particle: 𝐮=𝐔_𝑖+𝜔_𝑖×𝐱−𝐗_𝑖on𝜕𝑃_𝑖(𝑡). (10) Underlying assumption in this description is that the inertial and gravity forces inside the rigid ring may be safely neglected, which can be satisfied in this problem as well, in spite of the presence of the droplet impact. One way to check this applicability is to measure the particle Reynolds number and the drop Reynolds number Re_𝑑 and the particle Reynolds number Re_𝑝, are defined as Re_𝑑=𝜌_𝑑𝑉₀𝐷_𝑑𝜂_𝑑,Re_𝑝=𝜌_𝑑𝑉₀𝐷_𝑝𝜂_𝑑=Re_𝑑𝐷_𝑝𝐷_𝑑, (11) where 𝐷_𝑝 and 𝐷_𝑑 are the diameters of the particle and drop, respectively. In a typical inkjet printing process, the drop Reynolds number is less than 10 and, considering the drop and particle size ratio 𝐷_𝑝/𝐷_𝑑≪1, the particle Reynolds number becomes less than 𝑂(1), which facilitates the use of the rigid-ring description. Although not presented here, we tested another scheme, with which the rigid-body constraint is applied not only on the particle boundary but also in the particle interior, and the resultant rigid-body velocities were compared with results from the rigid-ring approximation. The two schemes were found to yield identical result.

The movement of the circular particle is given by kinematics 𝑑𝐗_𝑖𝑑𝑡=𝐔_𝑖,𝐗_𝑖∣_𝑡=0=𝐗_𝑖,0. (12) To determine the unknown particle rigid-body motion (𝐔_𝑖,𝜔_𝑖), one needs the balance conditions for the drag force and the torque on the particle boundary. In the absence of inertia and external forces or torques, particles are force-free and torque-free: 𝐅_𝑖=_{𝜕𝑃𝑖(𝑡)}𝑇𝝈⋅𝐧𝑑𝑠=0,_𝑖=_{𝜕𝑃𝑖(𝑡)}𝐱−𝐗_𝑖×(𝝈⋅𝐧)𝑑𝑠=0. (13) We did not use an artificial particle-particle collision scheme as in Glowinski et al. [14], as particle overlap and particle/wall collisions could be avoided through the use of a relatively small time step and sufficiently refined particle boundary discretization.

3. Numerical Methods

Jeong et al. [1] presented the weak form for capillary spreading of a particle-laden droplet using (i) the standard velocity-pressure formulation for the Stokes flow, (ii) the combined weak formulation of Glowinski et al. [14] for implicit treatment of freely suspended particles, and (iii) the level-set method with the CSS formulation [17] for the interface description between droplet and fluid. In the present work, an additional treatment is necessary to incorporate the inertia and convective terms in the Navier-Stokes equation for the description of impact and subsequent spreading behaviors. For this purpose, we employ the second-order Adams-Bashforth/Crank-Nicholson time stepping scheme to solve the Navier-Stokes equation [20]. In this scheme, the convective term is treated explicitly with the second-order Adams-Bashforth (AB2) method and the viscous term is treated implicitly with the Crank-Nicholson (CN) method. This time-stepping scheme for the Navier-Stokes equation can be written in a strong form as follows: 𝜌𝐮^𝑛+1−𝐮^𝑛3Δ𝑡=−𝜌2𝐮^𝑛⋅∇𝐮^𝑛−12𝐮^𝑛−1⋅∇𝐮^𝑛−1+12∇⋅𝝉^𝑛+1+∇⋅𝝉^𝑛−∇𝑝^𝑛+1+𝜌𝐠. (14) In the above equation, we introduce an extra stress tensor 𝝉 that includes the viscous stress 2𝜂^∗𝐃 and the interfacial stress 𝝉_𝑠.

The numerical method can be best interpreted with the time stepping procedure due to the AB2/CN time stepping for the momentum equation as below. At every time step, we take the following four step procedures.

Step 1. For given fluid velocity 𝐮^𝑛, we determine the level-set function 𝜑^𝑛+1 by solving the evolution equation of the level-set function: _Ω𝜓𝜑^𝑛+1−𝜑^𝑛Δ𝑡+𝐮^𝑛⋅∇𝜑^𝑛𝑑Ω−_𝑒Γin𝑒𝜓𝜑^𝑛−𝜑^𝑛_ext𝐮^𝑛⋅𝐧_𝑒𝑑𝑠=0. (15) To deal with the convection term, we employ the discontinuous Galerkin method in (15), in which 𝐧_𝑒 is the outward unit normal vector on the boundary of an element 𝑒, Γⁱⁿ_𝑒 is the part of the element boundary where 𝐮⋅𝐧_𝑒<0 (inflow), and 𝜑^𝑛_ext is the level-set function value in the neighboring upwind element. Although (15) is expressed in terms of the explicit Euler method, what is actually employed is the TVD-RK3 (total variance diminishing 3rd order Runge Kutta) method by applying this explicit Euler method three times for the single time step advance [17, 21].

Step 2. Given the particle velocity 𝐔^𝑛, the next particle positions 𝐗^𝑛+1 are updated (12).

Step 3. Given 𝜑^𝑛+1, the interfacial stress is evaluated 𝝉_𝑠^𝑛+1 from (6) using 𝛿(𝜑^𝑛+1).

Step 4. Given 𝜑^𝑛+1,𝝉_𝑠^𝑛+1, and 𝐗^𝑛+1, the weak form for the momentum and mass balances are solved with the constraint equation for the freely suspended particles. The weak form can be stated as follows: find (𝐮^𝑛+1,𝑝^𝑛+1,𝐔_𝑖^𝑛+1,𝜔_𝑖^𝑛+1,𝜆_𝑖^𝑛+1) such that −_Ω𝑝^𝑛+1()1∇⋅𝐯𝑑Ω+2_Ω2𝜂^∗𝐃𝐮^𝑛+1[𝐯]+1∶𝐃𝑑ΩΔ𝑡_Ω𝜌^∗𝐮^𝑛+1+⋅𝐯𝑑Ω_{𝑁𝑖=1𝜕𝑃𝑖}𝜆_𝑖^𝑛+1⋅𝐯^𝑛+1−𝐕^𝑛+1−𝝌^𝑛+1×𝐱−𝐗_𝑖^𝑛+1=𝑑𝑠_Γ𝑡𝐭_Γ⋅𝐯𝑑Ω−𝜌^∗_Ω32𝐮^𝑛⋅∇𝐮^𝑛−12𝐮^𝑛−1⋅∇𝐮^𝑛−1−1⋅𝐯𝑑Ω2_Ω2𝜂^∗𝐃𝐮^𝑛∶𝐃[𝐯]𝑑Ω+_Ω𝜌^∗()+1𝐠⋅𝐯𝑑ΩΔ𝑡_Ω𝜌^∗𝐮^𝑛1⋅𝐯𝑑Ω−2_Ω𝝉_𝑠^𝑛+1+𝝉^𝑛_𝑠∶𝐃[𝐯]𝑑Ω,_Ω𝑞∇⋅𝐮^𝑛+1𝑑Ω=0,_𝜕𝑃𝑖𝝁_𝑖⋅𝐮^𝑛+1−𝐔^𝑛+1−𝝎^𝑛+1×𝐱−𝐗_𝑖^𝑛+1(),𝑑𝑠=0𝑖=1,…,𝑁 (16) for all admissible functions (𝐯,𝑞,𝐕_𝑖,𝜒_𝑖,𝝁_𝑖)(𝑖=1,…,𝑁). In (16), the AB2/CN scheme is incorporated within the weak form and the rigid-ring constraint (10) has been combined as well by introducing the Lagrangian multiplier (𝑛=0) defined on the boundary of each particle. The Lagrangian multipliers can be (roughly) identified as traction forces acting on the particle by the fluid [13].

We have a few remarks on this four step procedures. At the first time step 𝜑⁰, the level-set function 𝐗⁰_𝑖 and particle positions (𝑖=1,…,𝑁)𝐮⁰ are set and the problem in Step 4 is solved to obtain velocities 𝐔⁰_𝑖 and 𝐮⁻¹ of particles and fluid. Then, at each time step, we repeat the four step procedures to advance to the next time step. The AB2/CN cannot be implemented at the first step as the fluid velocity at the previous time step 𝝉_𝑠 is not known and the explicit-Euler/CN has been used instead. As the level-set function evolves over time, the distance function property of (4) is no longer satisfied and the level-set function must therefore be reinitialized to restore the distance function property, because the accuracy of the distance function is crucial for the evaluation of the interfacial stress 𝐮. In this work, we reinitialize the level-set function by a direct geometrical redistancing method, following the work of Kim and Hwang [17]. We search for line segments along the interface on the subdivided element boundary and introduce a separate function, the distance function, whose function value is the shortest distance from any nodal position to the line segments along the interface. Details of the geometric construction of the distance function can be found in Kim and Hwang [17].

For the discretization of the weak form, we use a uniform regular quadrilateral element with the continuous biquadratic interpolation for the velocity 𝑝, a discontinuous linear interpolation for the pressure 𝜑, a discontinuous biquadratic interpolation for the level-set function 100×50, and the collocation point method the Lagrangian multipliers for the rigid-ring constraint. As the discretization method in the present study is the same as one of the author’s previous work, details of the implementation technique can be found in Hwang et al. [13], Kim and Hwang [17], Jeong et al. [1]. At every time step, we solve a large matrix equation from the discretization of the weak form (16) using a direct method based on a sparse multifrontal variant of Gaussian elimination (HSL/MA41).

4. Numerical Results

In this section, we present numerical results for the impact and spreading of a particulate droplet on a flat surface. A finite-element mesh of 𝑅^𝑖_𝑑=2[mm] is shown in Figure 4 along with the initial configuration of a droplet with 50 particles inside. The computational domain is of size [𝜇m], the initial drop radius is 𝜙 and the particle radius is 140 𝜙=0.245. The drop center is located initially at (5 mm, 1.9 mm) such that the bottom part of the drop is attached onto the solid surface from the beginning, as mentioned earlier. The particle fraction, denoted by 𝜌_𝑑, is defined as the ratio of the particle area to the total drop area and is fixed as 𝜂_𝑑 throughout this study. In the authors’ previous work [1], we investigated the particle-laden droplet spreading by the capillary force in the absence of inertia four different particle fractions (from zero to 0.954) and found that kinetic energy dissipation and local shear rate increase with the particle fraction, which leads to slow drop spreading even in the drop spreading processes governed by the capillary force. For this reason, we have focused in this work a relatively high particle fraction to clearly identify particle effects on drop deformation in this much more complex problem with effects of droplet impact with inertia.

Figure 4 

An example finite-element mesh with the initial drop interface. Particles, denoted as collocation points along the particle boundary, are distributed randomly inside the drop.

A particle is described by uniformly distributed collocation points on its boundary and in this problem ten collocation points are used to represent a particle. Material properties are roughly taken from the data of Nicolas [12]. The density 𝜌_𝑎 and viscosity 𝜂_𝑎 of the drop fluid are 400 kg⁄(cm³) and 0.01 Pa·s, respectively. As for the air region, we choose the density 𝜂_𝑑/𝜂_𝑎=100 and the viscosity 𝜌_𝑑/𝜌_𝑎=4000 such that 𝜍 and 60^∘. The interfacial tension coefficient 𝑔 is set to 2 mN⁄m, which is ten times smaller than the experimental data of Nicolas [12], due to the numerical stability problem and the equilibrium contact angle is set to 𝑔. The instantaneous artificial gravity 𝑉₀ ranges from zero to 4000, which correspond to 0 to 20.389 g with the gravitational constant 𝑔. The initial impact velocity Re_𝑑=𝜌_𝑑𝑉₀𝐷_𝑑/𝜂_𝑑 is the velocity at the droplet center, which is generated by the application of Ca=𝜂_𝑑𝑉₀/𝜁 for 0.0002 [sec]. In this case, the drop Reynolds number We=𝜌_𝑑𝑉²₀𝐷_𝑑/𝜁, the Capillary number Δ𝑡, and the Weber number 𝑔=5.869 are in ranges from 1.45 to 22.73, from 0.66 to 161.45 and from 0.45 to 7.1, respectively. The time step g used in this study is 0.0002 [sec] for (𝜙=0.245) g and 0.00005 [sec] for 20.389 𝐿.

Plotted in Figure 5 are consecutive droplet deformation patterns of a particulate droplet Re=5.869 under the capillary spreading (Figure 5(a)) and the spreading accompanied by the impact to illustrate overall behaviors of inertial spreading in comparison with the capillary spreading ( The corresponding movie files for this numerical experiment are accompanied with the paper as supplementary data: sup_movie_Figure  5(a).avi (capillary spreading) and sup_movie_Figure  5(b).avi (impact and spreading) available online at doi: 10.1155/2012/687961). (Note that the dimensions in all figures in this work are normalized by Ca=1.834). The dimensionless parameters in the impact problem are the drop Reynolds number We=10.765, the capillary number 𝑡=0.06, and the Weber number 𝑡=0.57. In the case of capillary spreading, as time passes, the drop spreads slowly and monotonically toward the equilibrium configuration. In contrast, inertial spreading looks more dynamic and indeed shows oscillation behavior: it reaches maximum spreading quickly at 𝑡=0.09 [sec] and then recoils back to near the original shape and the maximum recoil behavior is observed at 𝑡=0.51 [sec]. Subsequently, the drop spreads slowly and monotonically to the equilibrium state. To recognize the particle movement clearly, we colored three distinctive particles in black originally located near the bottom surface for both cases. As seen in Figure 5 as well as in the supplementary movies, the most evident motions of particles are generated by the impinging flow toward the surface in both cases, which is mainly driven by the spreading of the interface. In case of the impact spreading, the particle movement in the early fast spreading (before recoil) is found to be minor but it is observed apparently during the second spreading process after recoil. Particle redistribution is enhanced in the impact spreading due to high kinetic energy and active spreading in this case.

(a)

(b)

Figure 5 

The consecutive drop deformation shapes with particle distributions (a) for the capillary spreading and (b) the inertial impact spreading.

Next we investigate effects of particles on the droplet oscillation behaviors. Presented in Figure 6 are droplet interfaces at the instances of maximum spreading and maximum recoil for both the homogeneous fluid droplet and particulate droplet under the inertial impact spreading. Simulation conditions are taken the same as the previous inertial spreading of a particulate drop (Figure 5(b)). The homogeneous drop is observed to reach the maximum spreading at 𝑡=0.06 [sec] (Figure 6(a)) and the maximum recoil at 𝑡=0.6 s (Figure 6(b)), whereas the particulate drop reaches the maximum spreading slightly faster at 𝛽_𝑠 [sec] but the maximum recoil occurs later at 𝛽_ℎ than the homogeneous fluid drop case. In other words, the homogeneous drop spreads more and longer during the early spreading (before recoil) but also recoils faster than the particulate drop. This implies that the presence of particles reduces the mobility of the droplet spreading and recoil dynamics. Investigating the droplet interface at the maximum spreading and at the maximum recoil, one can observe that the homogeneous fluid drop spreads more and recoils more than the particulate drop, which indicates again the reduced mobility of the particulate drop.

(a)

(b)

Figure 6 

The droplet interfaces at the instances of (a) the maximum spreading and (b) the maximum recoils for a homogeneous fluid drop and the particulate fluid drops.

In order to provide physical insight of droplet spreading and oscillation, we present a high Reynolds number case (Re = 22.73, We = 161.45, and Ca = 7.1) as two additional supplementary movies: one for the pure fluid case (sup_movie_high_Re_purefluid.avi) and one with the presence of particles (solid fraction = 0.245, sup_movie_high_Re_suspension.avi).

For quantitative comparison of the droplet deformation, we present in Figure 7 the spread factor 𝛽_𝑠 and the splat height factor 𝛽_ℎ for three different cases discussed in the previous examples: (i) the capillary spreading of a homogenous fluid droplet; (ii) the inertial impact spreading of a homogenous fluid droplet; (iii) the inertial impact spreading of a particulate droplet. The spread factor is the droplet width normalized by the initial drop diameter and the splat height is defined as the droplet height divided by the initial droplet diameter. From Figure 7(a) for the spread factor, the homogenous droplet with impact is found to show the largest spread in the early inertial spreading phase, whereas the capillary spreading of the homogeneous drop does not show any oscillatory dynamics. The inertial spreading of the particulate droplet shows intermediate behaviors, which is consistent with the previous observation with Figure 6. The evolution of the splat height in Figure 7(b) shows monotonic decrease in case of capillary spreading of the homogeneous fluid droplet but fast drop followed by recoil back in case of the inertial impact and spreading. In the latter case, the amount of decrease in the splat height seems to be reciprocal to the particle fraction, which indicates again that particles restrain the oscillation of the drop. We remark that the evolution of the splat height in the homogeneous fluid droplet shows qualitative agreement with the data of Khatavkar [8], though not presented here.

(a)

(b)

Figure 7 

The comparison of time evolutions of (a) the spread factor 0.5(𝑢2+𝑣2) and (b) the splat height 𝛽ℎ for capillary and inertial spreading with the presence or absence of particles.

Suppressed oscillation due to the presence of particles may be related to energy dissipation, as discussed in Jeong et al. [1]. The presence of particles increases the shear rate and thereby the shear stress near the particles, which is expressed macroscopically by the increased bulk viscosity of the particulate fluid. The local energy dissipation can be expressed as the product of the increased bulk viscosity and square of “averaged” local shear rate (or the viscosity times square of increased shear rate due to presence of particles) and the observation of kinetic energy within a droplet during spreading may yield rough estimates for the amount of the energy dissipation. In this regard, we present the average kinetic energy in time in Figure 8, which is defined as 0.5(𝑢²+𝑣²) over the drop area including the interior of particles, for several Reynolds numbers with or without particles inside. (The particle velocity can be roughly identified by the average fluid velocity surrounding the particle.) One can observe completely distinct behavior of capillary spreading and the inertial spreading. In the capillary spreading, the kinetic energy increases due to the acceleration due to the capillary force, though slow, and then as the drop shape approaches the equilibrium configuration it vanishes to zero. On the other hand, there is a fast decrease of kinetic energy during the early spreading phase and then uprise of kinetic energy can be observed due to recoil of droplet. This means that a certain amount of kinetic energy is stored as the elastic energy of the interface during the first decrease in kinetic energy. Another observation is that the initial kinetic energy increases with the Reynolds number, as it should be. More importantly, one can see kinetic energy of the particulate drop decrease faster than that of the homogeneous fluid droplet in all cases, which supports our arguments of enhanced energy dissipation due to presence of particles, and the amount of uprise due to droplet recoil is again found to decrease with the particle presence. This increased energy dissipation seems to be responsible to the suppressed oscillatory behavior of the particulate droplet.

Figure 8 

The averaged kinetic energy in time during spreading for various flow conditions.

5. Conclusions

In this work, we presented a 2D numerical simulation technique for a particulate droplet spreading primarily driven by the inertial impact on a flat solid substrate and investigated particle effects on the particulate droplet dynamics such as spreading and recoil behaviors. Based on a finite-element method for the Navier-Stokes equation, the level-set method has been employed to describe the droplet/air interface dynamics and a DLM-like fictitious domain method has been introduced to describe freely suspended particles within the droplet. The combination of these techniques facilitates the use of a single Eulerian mesh in the entire computation.

We report that unlike the capillary spreading a particulate droplet reaches the maximum spreading quickly from the impact and recoils back to near original shape; and then spreads monotonically to an equilibrium state. During the course of impact and spreading, particle movements are found to be mainly driven by the latter monotonic spreading. In addition, we observed that the presence of particles reduces the mobility of the droplet spreading and recoil dynamics that invokes the reduced amount of the spread factor at the maximum spreading and slower deformation during spreading and recoil. Investigating changes in the average kinetic energy in time, we hypothesized that the increased energy dissipation due to the presence of particles would be responsible for the reduced mobility of the particulate droplet.

Acknowledgments

This work was supported by the National Research Foundation (NRF) of Korea funded by the Ministry of Education, Science, and Technology (2011-0031383 and 2010-0021614).

Supplementary Materials

One for the pure fluid case (sup movie high Re purefluid.avi) and one with the presence of particles (solid fraction = 0.245, sup movie high Re suspension.avi)

1. Supplementary video 1
2. Supplementary video 2
3. Supplementary video 3
4. Supplementary video 4