Forced oscillations in sessile droplets can be exploited in electrowetting mixing of fluid fractions. The necessary complex flows and large shape deformations require a numerical investigation of fluid dynamics in the transient regime. We provide a means to characterize oscillations qualitatively and quantitatively with the goal to examine and to classify flow patterns occurring inside. A superposition of different harmonic excitation patterns gives the possibility to control the convective flow. In this investigation, we apply a generic and accurate multi-phase smoothed particle hydrodynamics model to a two-dimensional three-phase flow and consider an oscillating droplet sitting on a substrate and immersed in a fluidic phase. These vibrations are investigated in two ways: the analysis of a step response due to an abrupt change of the contact angle is applied to identify the resonance frequencies. Secondly, the time evolution of the shape of the droplet in terms of harmonic functions is determined. Their amplitudes are examined in the time and frequency domain. This gives the possibility to relate resonance frequencies to mode shapes and to detect a coupling between them. Our approach is successfully applied to different numerical case studies.
Oscillations in liquid, viscous droplets have been considered both analytically and numerically [1–4]. The corresponding solution of an eigenvalue problem was limited to small shape deformations of the linearized Navier–Stokes equation. With the prospect of lab-on-a-chip applications, forced droplet oscillations by the electrowetting (EW) effect have been studied experimentally [5–7]. The droplets perform self-excited oscillations that were compared for small amplitudes with the aforementioned linearized models. An investigation on droplet shapes with respect to EW actuation and surface tension models can be found in Lienemann et al.  where suggestions are given on how to improve EW-induced actuation.
For an accurate analysis of the transient regime, including large amplitudes and different transport properties, a multi-phase capable numerical model is needed. For simulating sessile droplets, an adjustable surface tension between the different phases allows for setting an apparent contact angle of the liquid phase.
An investigation of droplet oscillations requires stability over a large number of time steps, computational efficiency, as well as long-term numerical stability of the interface. Thus, conventional Eulerian methods like the finite-element method or the finite-volume method are challenged: rapid changes in the topology require frequent and expensive remeshing. A further transport equation for the phases like in a phase-field model has to be solved and smears the interfaces in the long run: phase fractions are transported cell-wise and spread by numerical diffusion. Lagrangian, meshless methods, on the other hand, are well suited for the requirements above. The discretization points, the particles, carry the information to which phase they belong and allow for sharp interfaces, even in the long run. From the different meshless methods, smoothed particle hydrodynamics (SPH) seems to be especially attractive for being relatively well explored and cheap.
Oscillations of a two-dimensional droplet with a given contact angle are excited by abruptly changing the wetting properties of the substrate. Thereby, all eigenfrequencies are excited, the eigenspectrum is obtained and an estimate of the resonance frequencies can be given. For measuring mode shapes, the surface of the droplet is represented in terms of spherical harmonics. We calculate the mode in the time and frequency domain, relate the resonance frequencies to mode shapes and identify mode coupling.
2. Model and analysis
We compute a multi-phase flow of a liquid and a vapour phase on a solid substrate by solving the Navier–Stokes equations with a multi-phase SPH model. By introducing different surface energies between the phases involved, corresponding contact angles arise. This is illustrated in figure 1 for contact angles of 60° and 120°. By changing the surface energy between solid and liquid phase, vibrations are induced in the droplet.
(a) Governing equations
We consider the isothermal and incompressible Navier–Stokes equations in Lagrangian form 2.1and 2.2with velocity v, density ρ, pressure force F(p), surface tension force F(sf), viscous force F(v) and gravitational acceleration g. For incompressible media, the viscous part of the momentum equation with a dynamic viscosity η simplifies to 2.3
Incompressibility is approximated by a stiff equation of state (EOS), which relates pressure p to density and allows us to keep density fluctuations small , 2.4where p0 represents a reference pressure, ρ0 a reference density and γ depends on the phase, where a usual choice is 7 for liquids and 1.4 for gases. The continuum surface force (CSF) model from Brackbill et al.  gives 2.5with the surface tension coefficient α, the curvature κ and the surface normal n. The surface delta function δsf describes the narrow transition band between two phases and in the limit approaches the Dirac delta function representing a spatial phase jump. The equilibrium contact angle of the liquid phase that arises is given by the Young equation, which relates surface energies to a contact angle, 2.6where γij is the surface energy between two phases i and j, ϕ is a contact angle and s, l, v, respectively, label the solid, liquid or vapour phase.
(b) Smoothed particle hydrodynamics multi-phase model
For zeroth-order consistency, we use a normalized smoothing function , 2.7
where W is a generic SPH smoothing kernel, h is the smoothing length and is a global sum contributing only within the support of W. The usual SPH density equation  then reads 2.8so that the density of a particle i does not depend on the mass mj of surrounding particles, but only on its specific volume . For continuity at interfaces, we use an expression suggested by Hu & Adams  for a particle-averaged spatial derivative 2.9of a quantity Φi, where , is an interparticle average and eij is a unit vector pointing from particle i to particle j.
The mass in equations (2.1) and (2.8) is inherently conserved by the Lagrangian nature of our approach. The pressure force with the interparticle pressure is discretized as 2.10Secondly, the shear stress contributes to the momentum equation and for different viscosities ηk and ηl in different phases k and l, the interparticle-averaged shear stress can be written as 2.11when assuming that the phase interface is centred between two particles. For incompressible flow, the following expression is found : 2.12
Furthermore, the surface tension force of the CSF model contributes to the momentum equation. A colour index is defined for each phase. A corresponding colour function defined on each particle pi is 1 if particle i belongs to phase s, otherwise it is 0. The gradient of the colour function for a particle i in phase k between two phases k and l can be written as 2.13where an interparticle average is avoided. Considering equation (2.5), a straightforward SPH discretization for the calculation of the interface curvature κ=∇⋅((∇C)/(|C|)) is cumbersome for high curvatures  and a tensorial description that avoids direct computation of second derivatives is used instead. Hu & Adams  suggest a multi-phase compatible description of surface stress between two phases k and l 2.14where I is the unit tensor, d the dimensionality, is the global surface stress and F(sf)=−ρ−1∇Π(sf), which we adopt.
(c) Mode shape detection and modal analysis
We express the radius of a droplet with respect to its mass centre, which is projected onto the substrate’s surface, by a superposition of spherical harmonics functions 2.15with the inclination ϑ, the azimuth φ, a spherical harmonics function Y lm and its attributed amplitude clm. For our analysis, we interprete two-dimensional functions as φ-symmetric. A particle is represented in spherical coordinates with the components r,ϑ and φ. The set of particles forming the droplet is situated on a half sphere with the angular range . We divide the sphere into equally sized bins ϑi and φi, i=1,…,n, with the corresponding particle subsets and . We define a radius function to approximate the distance of a sphere segment of the droplet’s surface from its centre 2.16where Ξr returns the radial component of a particle p. The coefficients clm of the individual spherical harmonics for the projection of equation (2.16) and expansion with equation (2.15) are calculated by 2.17with being the conjugate complex of Y lm.
To determine the eigenspectrum, the initial surface energy αsl,i is changed abruptly to αsl,s corresponding to a step excitation and the step response is measured. The approach of applying the Fourier transformation on the trajectory ri(t) of individual particles was useful for obtaining the eigenspectrum when integration errors were substantially lowered by tremendously reducing the time-step length. Transforming the expansion of the surface of the droplet by equation (2.15) is numerically more efficient. Only the noise inherent in the trajectories of surface particles contributes to errors and as a further benefit, the number of expensive Fourier transformations is thereby decreased.
3. Numerical results
The simulations in this section were performed on a rectangular periodic domain [−0.5;0.5]×[−0.4;0.4] filled with 2880 particles having initial positions on a regular lattice, a quintic smoothing kernel with the smoothing length and a support of 3h, equal viscosities η=0.05 and densities ρ0=1 for all phases and an artificial speed of sound of vs=100. Initial surface tension coefficients were chosen as αsl=αlv=αsv=1. All simulations ran with a time step of Δt=10−6 and for about 4×106 time steps. For all simulation runs, we calculated a maximum Reynolds number of and Weber numbers in the range of We=[0.5;1.5].
(a) Step response
We change the solid–liquid surface energy abruptly from 0.5 to 1.5 which corresponds to steady-state contact angles of 60° and 120°, respectively. The surface of the droplet is described by equation (2.16). The expanded droplet’s radius from equation (2.15) is then Fourier transformed.
In figure 3a, we observe some peaks in the spectrum for lower frequencies at about f≈2 Hz,f≈5 Hz and f≈7 Hz. In the calculated sample, numerical noise appears to inhibit the detection of lower amplitude eigenfrequencies. A higher sampling rate, as well as repeated excitation, might improve results. The latter is presented in figure 3b by a square wave excitation with identical system configuration and improved methodology.
(b) Mode shapes
Now, we consider individual mode shapes in the time and frequency domain. The droplet is excited with different patterns, where αsl,0 corresponds to a constant surface energy term, αsl,A to an amplitude, and sgn stands for the signum function, 3.1and 3.2All steady surface energies were set to αsl,0=αlv=αsv=1, corresponding to a contact angle of 90°. The excitation frequency was set to f=5 Hz.
Owing to enforced quasi-incompressibility, the Y 00 mode is approximately constant, presented for both patterns in figure 3b: the small deviations reflect the stiffness of the EOS. For symmetry reasons, odd modes do not change substantially. Figure 4a shows first interesting amplitudes c20,c40 and c60 for a harmonic excitation by equation (3.1), and figure 4c for the square wave excitation by equation (3.2). In both figures, the c20 amplitude slowly increases, indicating resonance behaviour. This corresponds to the peaks detected in §3a and observed in the eigenspectrum in figure 3a.
Figure 4b,d shows the clm in the frequency domain. The peak at f=5 Hz for all modes corresponds to the excitation frequency. For both excitation patterns in the clm spectra, a contribution to resonance is visible for c20 at f≈2 Hz. The amplitude for the square wave excitation is however slightly higher, which might be caused by the higher amount of energy brought into the system. A further correspondence to the eigenspectrum is observable at a frequency of f≈7 Hz, where different modes show contributions.
We finally relate our results to an analytical model , which gives resonance frequencies for small amplitude oscillations of cylindrical droplets in an external fluid. ρd and ρe are corresponding densities, l stands for the mode number, σ for the droplet’s surface tension and r0 for the radius of the unperturbed sphere. For the meaningful second mode, the model predicts f=10.54, which is in fair agreement with our result from a damped system.
4. Concluding remarks
We have presented an SPH-based normal mode analysis of a sessile, oscillating droplet and provided a means to extract and to represent its shape in terms of spherical harmonics functions. The presented approach for estimating resonance frequencies and measuring mode shapes in the time and frequency domain is broadly applicable, as the multi-phase scheme applied allows for different transport properties in the individual phases. Numerical case studies were performed and compared with resonance frequencies in order to relate individual mode shapes. An extension to complex flows is straightforward and makes possible meaningful parameter studies for industrial applications beyond analytical scope.
The authors acknowledge complete funding by the DFG via the project Electrowetting-Simulation mit Partikelmethoden (grant no. LI 1831/1-1).
One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.
- This journal is © 2011 The Royal Society