1. Introduction
Dynamic wetting of a solid substrate by a liquid film is a common occurrence. In natural systems, it appears in surfactant films that wet the lungs to promote breathing and the eyes to keep them moist (Holly & Lemp Reference Holly and Lemp1971; Veldhuizen & Haagsman Reference Veldhuizen and Haagsman2000), or in the spreading of water drops on solid surfaces, such as raindrops on windows (Dussan V. & Davis Reference Dussan V. and Davis1974; Dussan V. & Chow Reference Dussan and Chow1983). In technological systems it is harnessed for actuating microfluidic platforms (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Atencia & Beebe Reference Atencia and Beebe2005; Whitesides Reference Whitesides2006), for cooling electronic circuits (Amon et al. Reference Amon, Murthy, Yao, Narumanchi, Wu and Hsieh2001; Bar-Cohen, Arik & Ohadi Reference Bar-Cohen, Arik and Ohadi2006), for desalination (Fletcher, Sernas & Galowin Reference Fletcher, Sernas and Galowin1974) and for a variety of manufacturing applications (Fendler Reference Fendler1996; Nagayama Reference Nagayama1996; Wang et al. Reference Wang, Zheng, Li, Huck and Sirringhaus2004). Many different mechanisms govern dynamic wetting, where capillary forces at the free surface of a liquid film, or at a meniscus, contribute appreciably to the dynamics of the film. In addition to capillary forces, dynamic wetting may be powered by gravity, thermal or solutal Marangoni effects or electrokinetics, among other driving mechanisms; see Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997) or Craster & Matar (Reference Craster and Matar2009) for comprehensive reviews. In the present work, we focus on the contributions of mechanical vibrations in the solid substrate, particularly MHz-frequency Rayleigh surface acoustic waves (SAWs), to dynamic wetting (Morozov & Manor Reference Morozov and Manor2018). We study particularly the effects of bulk acoustic streaming – mass transport in the bulk fluid that results from variations in acoustic intensity therein. While extensively observed and described in the scientific literature (Brunet et al. Reference Brunet, Eggers and Deegan2007, Reference Brunet, Eggers and Deegan2009; Yeo & Friend Reference Yeo and Friend2014), we are not aware of a model describing the dynamics of such a system that further accounts for the capillary and gravitational stresses that shape the liquid/air interface for the case when drops are sufficiently thick (millimetre scale, therefore much thicker than the SAW wavelength). Here, we use thin film (long wave) theory alongside physical experiments to develop a better understanding of the dynamics of this problem.
There are several different mechanisms by which acoustic waves or other mechanical waves, at a solid boundary or at the free surface of a liquid, generate streaming – steady flow along the wave’s path. A mechanical wave at a boundary invokes periodic viscous flow within the viscous penetration length
$\delta \equiv \sqrt {2\mu /(\rho \omega )}$
away from the interface, where
$\mu ,\,\rho$
and
$\omega$
are the viscosity and density of the liquid, and the wave angular frequency, respectively. Examples in which this problem was considered include the work by Rayleigh on a standing acoustic wave grazing a solid bed (Rayleigh Reference Rayleigh1884); by Schlichting (Reference Schlichting1932) on standing vibration of waves of infinite and finite (respectively) wavelength in a solid along its surface; by Longuet-Higgins (Reference Longuet-Higgins1953) on shallow ocean waves; and by Manor, Yeo & Friend (Reference Manor, Yeo and Friend2012) on propagating Rayleigh (surface acoustic) waves in a solid substrate. Convective contributions due to surface waves invoke a drift of liquid mass, which does not attenuate away from the solid surface. This is known as Schlichting streaming (Yeo & Friend Reference Yeo and Friend2014), and its component far from the boundary was further coined the ‘Rayleigh law of streaming’ by Lighthill (Reference Lighthill1978) in recognition of Rayleigh’s work on the drift resulting from the presence of standing acoustic waves in fluid near a solid boundary.
A separate mechanism, the Eckart streaming, appears regardless of the presence of a boundary. This is an acoustic streaming that occurs in the bulk fluid, and results from variations in the intensity of sound or ultrasound waves in the fluid. In the classic work by Eckart and in many subsequent studies, acoustic streaming is ascribed to the viscous attenuation of travelling waves in the fluid (Eckart Reference Eckart1948; Nyborg Reference Nyborg1952; Lighthill Reference Lighthill1978). This attenuation results in spatial variations of the wave intensity in the fluid, and in corresponding spatial variations in convective Reynolds stresses. The consequence is a secondary flow, whose steady component at long times is the acoustic streaming, characterised by an intense vortical flow field that appears in the bulk fluid. This is the main acoustic mechanism for the actuation of fluid in micro-channels (Wixforth Reference Wixforth2003; Wixforth et al. Reference Wixforth, Strobl, Gauer, Toegl, Scriba and von Guttenberg2004; Yeo & Friend Reference Yeo and Friend2014; Fasano et al. Reference Fasano, Diez, Manor, Kondic and Cummings2025) and in drop microfluidics (Guttenberg et al. Reference Guttenberg, Rathgeber, Keller, Rädler, Wixforth, Kostur, Schindler and Talkner2004; Brunet et al. Reference Brunet, Eggers and Deegan2007, Reference Brunet, Eggers and Deegan2009, Reference Brunet, Baudoin, Matar and Zoueshtiagh2010). It usually dominates contributions to the flow in the bulk fluid arising from the boundary layer-type acoustic streaming mechanisms, such as Schlichting or Rayleigh streaming described above. An additional effect related to the interaction of acoustic waves with a surface, the acoustic radiation pressure (Hertz & Mende Reference Hertz and Mende1939; Hamilton & Blackstock Reference Hamilton and Blackstock1998), may be also relevant in the presence of interfaces, such as the fluid–air surface. This effect has been discussed extensively in the literature (King Reference King1934; Campbell Reference Campbell1970; Chu & Apfel Reference Chu and Apfel1982; Shiokawa, Matsui & Ueda Reference Shiokawa, Matsui and Ueda1989; Hasegawa et al. Reference Hasegawa, Kido, Iizuka and Matsuoka2000; Borgnis Reference Borgnis2010; Karlsen, Augustsson & Bruus Reference Karlsen, Augustsson and Bruus2016; Rajendran, Aravind & Subramani Reference Rajendran, Aravind and Subramani2023), and shown capable of deforming and displacing soft interfaces and of governing the dynamics of films (Biwersi, Manceau & Bastien Reference Biwersi, Manceau and Bastien2000; Alzuaga, Manceau & Bastien Reference Alzuaga, Manceau and Bastien2005; Issenmann et al. Reference Issenmann, Wunenburger, Manneville and Delville2006; Rajendran et al. Reference Rajendran, Jayakumar, Azharudeen and Subramani2022; Marcos et al. Reference Marcos, Li, Fasano, Diez, Cummings, Manor and Kondic2025). Altshuler & Manor (Reference Altshuler and Manor2015) consider the problem of a partially wetting water drop atop a MHz frequency propagating vibration in the substrate, in a set-up similar to the one considered here, although the reduced viscosity and increased surface tension lead to reduced viscous dissipation of the leaky SAW in the fluid, and a necessity to include disjoining pressure to model liquid–solid interactions. Horesh et al. (Reference Horesh, Khaikin, Karnilaw, Zigelman and Manor2019) further consider a vertical set-up where gravity becomes important and show that acoustic radiation pressure is negligible for the partially wetting water due to the curvature of the meniscus, which does not support acoustical resonance.
While the literature discussing the interaction of acoustic waves with fluids is clearly extensive, in what follows, we discuss briefly just a few works that are particularly relevant in the present context, in which a SAW travelling in a solid substrate drives flow in a neighbouring fluid. In 1970, Campbell & Jones (Reference Campbell and Jones1970) simulated a SAW in a piezoelectric solid in contact with an ideal (inviscid) fluid. Then, in the late 1980s, Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989) used the ansatz of a harmonic acoustic waveform in the solid based on ideas developed by Nyborg (Reference Nyborg1952), converted the approach adopted by Campbell & Jones (Reference Campbell and Jones1970) to analytic expressions and calculated the resulting acoustic forcing in the adjacent inviscid fluid due to attenuation of SAW in the solid. More recently, Vanneste & Bühler (Reference Vanneste and Bühler2011) revisited the work of Campbell & Jones (Reference Campbell and Jones1970) in the context of viscous fluids and calculated numerically the corresponding acoustic forcing in the fluid.
In the present paper, we focus on developing a model that can be used to formulate a long-wave theory applicable to the free surface evolution of thin viscous films and drops due to the presence of a SAW propagating in an underlying solid substrate. To facilitate such a development, we simplify the previously considered approaches in order to isolate the dominant physical effects governing the dynamics of the film. In particular, we simplify the approach used by Campbell & Jones (Reference Campbell and Jones1970) by assuming damped harmonic waves in the solid. Furthermore, we compare the previous approaches by Vanneste & Bühler (Reference Vanneste and Bühler2011) and Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989) to quantify the importance of viscous effects, which were included by the former authors but ignored by the latter. We follow earlier attempts (Brunet et al. Reference Brunet, Baudoin, Matar and Zoueshtiagh2010) and focus on the contribution of acoustic streaming to the drop dynamics. In principle, one could also consider including the effects of the acoustic radiation pressure at the free surface, as was done in the recent work of Marcos et al. (Reference Marcos, Li, Fasano, Diez, Cummings, Manor and Kondic2025), who used a Monte Carlo-type approach to study the dynamics of oil–water films under acoustic forcing. However, in the present case, the main contribution to drop displacement is the acoustic streaming, which efficiently transfers power from the SAW to the translational motion of the drop; hence, we do not include acoustic radiation pressure here. This point will be discussed in further detail in § 3.
The rest of this paper is structured as follows. As a motivation, in § 2 we discuss simple physical experiments that demonstrate the well-known phenomenon of thick oil films and drops displacing under the action of MHz-frequency Rayleigh-type SAWs. These experiments illustrate the process of dynamic wetting and drop displacement powered by bulk acoustic streaming. We then derive the theoretical model that can explain such observations in § 3. Our theoretical results are presented in § 4: we discuss computational results in § 4.1, before analysing a simplified travelling–wave-type model in § 4.2. Sections 3 and 4, together with the appendices, constitute the main novel contributions of our work. Section 5 is devoted to the summary and conclusions. The appendices include (Appendix A) details of the theoretical results, (Appendix B) various simplifications and approximations including the inviscid limit and (Appendix C) the outline of the numerical approach implemented. Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10802 provide videos of selected experiments and animations of representative computational results.
2. Experiment
The experiments described here are illustrative and portray the well-known phenomenon in which MHz-frequency SAWs are used to displace drops and thick films of liquid along the path of the SAW. The purpose of these experiments is to demonstrate the system under investigation, highlight key features and trends and provide results that will be used to motivate the theory and simulations presented in the following sections.
Figure 1 shows the experimental set-up (A) and a typical experimental result (a–e). We generate a propagating 20 MHz frequency SAW by applying a same-frequency sinusoidal voltage signal to a piezoelectric actuator – a SAW device. The actuator comprises a 5 nm titanium/1
$\unicode{x03BC}$
m aluminium interdigitated transducer (IDT; from which the SAW emanates) fabricated atop lithium niobate (LiNbO
$_3$
, Roditi International, UK) by standard lift-off photolithography. The substrate used for the SAW device is 11 mm
$\times$
24 mm in size, 0.5 mm thick, 128
$^\circ$
Y-cut, X-propagating, single-crystal piezoelectric lithium niobate, where X and Y are crystal axes (Campbell & Jones Reference Campbell and Jones1968). The actuator is integrated into the external electrical signal using pogo pins (BC201403AD, Interconnect Devices, Inc.) assembled in a 3D-printed elastomeric stage, which holds the actuator and is connected to a signal generator (R&S SMB100A microwave signal generator) and amplifier (model A10160, Tabor Electronics Ltd.). Prior to the experiment we clean the SAW device using four different solvents: acetone (AR-b, 99.8 %, 67-641, Bio-Lab Ltd.), 2-propanol (AR-b, 99.8 %, 67-63-0, Bio-Lab Ltd.), ethanol (CP-p, 96 %, 64-17-5, Bio-Lab Ltd.) and water (HPLC plus, 7732-18-5, Sigma-Aldrich), before drying the actuator using compressed air. We place a paper cylinder soaked with glycerol at the far end of the actuator from the IDT to absorb the SAW and prevent reflections.
(A) Upper schematic view of the experimental set-up: SAW propagates from the IDT until it reaches the acoustic absorber (comprising glycerol-soaked paper placed on the actuator under and to the right of the needle). The needle is of known diameter (510
$\unicode{x03BC}$
m), placed for identifying spatial resolution in the images. (a–e) Successive snapshots are taken from an experiment monitoring the flow of a silicone oil film. During the experiment, a drop of silicone oil is placed on the horizontal surface (a), it deforms due to the application of SAW (b) and moves in the direction of SAW propagation (c–e); time is shown in seconds, and the vertical lines serve as a reference.
In our experiments, we represent the SAW strength (intensity) by measuring the normal displacement amplitude,
$A_{{n}}$
, at the solid surface before any attenuation in the liquid, so that
$A_{{n}}$
is the amplitude of the transverse component of the SAW in the substrate. We measure
$A_{{n}}$
over a surface of
$1\times 1\,\text{mm}^2$
approximately
$2$
mm away from the SAW actuator using a scanning laser Doppler vibrometer (MSA-500, Polytech). We control this amplitude by changing the applied voltage,
$V$
, of the signal generator at the surface of the actuator: consistent with the literature (Ballantine et al. Reference Ballantine, White, Martin, Ricco, Zellers, Frye and Wohltjen1996), there is a linear relationship between
$A_{{n}}$
and
$V$
; see Li et al. (Reference Li, Marcos, Fasano, Diez, Cummings, Kondic and Manor2025) for more details.
We place an
$8\,\unicode{x03BC}$
l (8 mm
$^3$
) drop of silicone oil (50 cSt, 378 356, Sigma-Aldrich) atop the actuator using a pipette, approximately 3.5 mm away from the IDT, and introduce electrical signal at different voltage levels to induce motion of the oil film, which dynamically wets the solid substrate along the path of the SAW. Figure 2(a–c) shows a time lapse (side view) of the typical silicone oil dynamics, captured using the framework and camera of a goniometer (Data Physics; OCA 15Pro). Figure 2(d) shows a top view obtained using a camera (EOS R5, Canon) with a macro lens (RF 100 mm F2.8L MACRO IS USM, Canon). To confirm that the observed dynamics is indeed induced by SAW, we also carried out control experiments, in which SAW is absent and all other parameters are kept the same. The oil drop in this experiment, which was repeated three times, slowly spreads concentrically, as expected in the absence of a directional forcing mechanism (figure not shown for brevity).
Two sets of images from experiments monitoring silicone oil drops powered by SAWs of three different measured amplitudes
$A_{{n}}$
(1.2–1.4 nm), where (A–C) silicone oil starts from the same location, moves in the direction of the SAW and reaches different distances (a–c) after 4.09 s. The corresponding side-view videos at different acoustic power levels are used to analyse the speed and profile of silicone oil during the movement. (D, d) Show top-view snapshots (corresponding to (C, c)) of the oil. A thin film of oil, hardly visible in the side view, can be seen more clearly behind the main body of silicone oil in the top views.
Time evolution of (a) the front position and (b) the maximum drop height for several values of
$A_{{n}}$
, for silicone oil drops of volume
${\mathcal{V}}_{{d}}=8\, \unicode{x03BC}$
l and kinematic viscosity
$\nu =50$
cSt. The symbols correspond to the measurements, and the lines are simple fits to guide the eye. We expect errors of
$\pm 50\,\unicode{x03BC}$
m in the values of the film height due to the limited resolution of the side-view image and tracker software. Note that
$t=0$
corresponds to the time instant at which SAW is applied, which is 2 s after deposition of the oil. A video showing the evolution at various SAW amplitudes is available as a supplementary material (movie 1).
In the experiments with SAW actuation we employ the side-view video, recorded by a goniometer camera, to capture the position of the advancing three-phase contact line of the SAW-actuated oil film using the open-source software Tracker (Brown & Cox Reference Brown and Cox2009) (see figure 3
a for typical data). We also measure the maximum height of the moving oil film using Tracker (typical results in figure 3
b). Figure 3(a) shows that, as expected, the larger the SAW amplitude
$A_{{n}}$
, the faster the oil film moves due to the stronger induced leakage wave (Shiokawa et al. Reference Shiokawa, Matsui and Ueda1989). Another notable aspect of SAW-induced dynamics is the contact line speed for early times (
${\lt } 1$
s): initially, it is rather small, but later increases and then remains constant. This change of contact line speed appears to be a consequence of the initial change of shape of the drop (see also figure 1); while this is taking place, contact line speed is small. After the initial deformation, the drop shape remains unchanged, translating approximately uniformly. Figure 3(b) shows consistent results for the maximum drop height, with an increase in height for early times, and constant values for longer times. Note that this uniform translation and constant film height differ significantly from the observed behaviour for spreading under a body force such as gravity. These experimental observations will be explained and discussed within the framework of our theoretical model, which is described in the next section. We note that in our experiments, we do not observe drop oscillations as in the experiments using water drops (Brunet et al. Reference Brunet, Baudoin, Matar and Zoueshtiagh2010), possibly due to the much higher viscosity of the silicone oil.
In the experiments measuring the asymptotic front speed, we place a drop of silicone oil of various volumes (
$4/8/16\,\unicode{x03BC}$
l) and viscosities (50 cSt, 378 356, Sigma-Aldrich; 100 cSt, 378 364, Sigma-Aldrich; 500 cSt, 378 380, Sigma-Aldrich) on the SAW actuator in a similar manner. We then apply different levels of SAW to induce the movement of the oil film. The initial frames of recorded side-view videos are processed by public domain image processing software ImageJ to mark a known distance of 5 mm (see the red lines in figure 2(A–C) and movie 1, movie 2, movie 3 in the supplementary materials). The times at which the drop fronts pass these 5 mm markers are recorded. Since the marked area is positioned rather far from the IDT, we assume the drops have by this time settled to uniform translation, and we use these markers to calculate the average speed of the drops (of various volumes/viscosities); see figure 4(a–b). Figure 4(a) shows that front speed decreases as viscosity increases, assuming a fixed amplitude of the SAW. Interestingly, this decrease in front speed does not scale linearly with viscosity (as is the case for drops under the influence of only gravity and capillarity, i.e. drops sliding down a ramp or spreading on a flat substrate under capillary forces) due to the nonlinear dependence of the acoustic terms on the viscosity, which is discussed in the model formulation in § 3. Figure 4(b) shows that the front speed increases as drop volume increases, an effect that is magnified for larger SAW amplitudes.
Asymptotic front speed of drop profiles for several values of the SAW amplitude
$A_{{n}}$
obtained using (a) silicone oil drops of volume
${\mathcal{V}}_{{d}}=8\,\unicode{x03BC}$
l and different kinematic viscosities
$\nu$
; and (b) different volumes
${\mathcal{V}}_{{d}}$
of silicone oil of the same viscosity
$\nu =50$
cSt. We present two example movies comparing typical evolution at various kinematic viscosities and volumes in the supplementary material (movie 2, movie 3).
3. Model
In our model development, we follow the approach of Vanneste & Bühler (Reference Vanneste and Bühler2011) and Campbell & Jones (Reference Campbell and Jones1970), implementing simplifications appropriate to our experiments. One important point is that, for the considered parameters, the dominant contributions to the Reynolds stress (responsible for the force driving the flow) arise from both the SAW attenuation in the solid, and from viscous damping of the leaked wave. The latter effects were not included by Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989); we will comment in particular in Appendix B on the consequences of this difference.
Schematic of a liquid drop (blue curve) driven by a leaky SAW (LSAW) of amplitude
$A_{{n}}$
. The SAW travels from left to right in the solid substrate. The drop extends from the rear to the front contact lines at
$x_{{r}}(t)$
and
$x_{{f}}(t)$
, respectively. The thin red line represents the amplitude
$A_{{n}}$
of the SAW at the substrate and shows how it is attenuated due to the presence of the liquid (the black dotted line is the envelope of the decaying amplitude in the liquid). The Rayleigh angle is denoted by
$\theta _R$
. The thickness
$h^\ast$
and the positions
$x^{\ast }_{1,2}$
are defined later in the text, see (4.3).
We consider a two-dimensional Newtonian fluid layer occupying a domain in the
$(x,z)$
-plane, traversed by a sound wave generated by a Rayleigh-type SAW travelling in a neighbouring solid; see figure 5 for a sketch. We assume that the associated acoustic field Mach number,
$M_a \equiv \omega A_{{n}}/c$
(where
$\omega$
,
$A_{{n}}$
and
$c$
are the SAW angular frequency, characteristic normal displacement amplitude at the solid surface and the phase velocity of the acoustic field in the fluid, respectively) is small (
$M_a \ll 1$
), and that the velocity (
$\boldsymbol{v}$
), pressure (
$p$
), and density (
$\rho$
) fields in the fluid can be written as
\begin{align} \boldsymbol{v} &= {\boldsymbol{v}}_0+{\boldsymbol{v}}_1 + {\boldsymbol{v}}_2+\cdots = {\boldsymbol{v}}_0+M_a \, \bar {\boldsymbol{v}}_1 + M_a^2 \, \bar {\boldsymbol{v}}_2+\cdots , \nonumber \\ p &= p_0 + p_1 + p_2+\cdots = p_0 + M_a \, \bar {p}_1 + M_a^2 \, \bar {p}_2+\cdots , \\ \nonumber \rho &= \rho _0 + \rho _1 +\rho _2+\cdots = \rho _0 + M_a \, \bar {\rho }_1 +M_a^2 \, \bar {\rho }_2+\cdots , \end{align}
where the leading-order components of these fields are associated with quiescent liquid (
${\boldsymbol{v}}_0=0$
,
$p_0$
constant) of constant ambient density
$\rho _0$
; the first-order corrections (
${\boldsymbol{v}}_1,\,p_1,\,\rho _1$
), represent an oscillatory flow associated with the acoustic field generated by the SAW that average to zero over the period of the acoustic forcing; and the second-order corrections (
$\boldsymbol{v}_2,\,p_2,\,\rho _2$
) represent the induced steady streaming flow. The assumption implicit in the above expansions in the Mach number is that the barred quantities are comparable in size to the corresponding leading-order contributions. The problem is governed by the Navier–Stokes and continuity equations,
\begin{align} &{\partial (\rho \boldsymbol{v})}/{\partial t} + \rho \left ( \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{v} + \boldsymbol{v} \,\boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \rho \boldsymbol{v} \right )\nonumber\\&\quad = -\boldsymbol{\nabla }p + \mu {\nabla} ^2 \boldsymbol{v} + \left ( \mu _b +{\mu }/{3} \right ) \boldsymbol{\nabla }\left ( \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v} \right ) + \rho \boldsymbol{g}, \end{align}
\begin{align} &\qquad\qquad\quad {\partial \rho }/{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\rho \boldsymbol{v}\right )=0, \end{align}
where
$t$
is time;
$\mu$
,
$\mu _b$
are the shear and dilatational coefficients of the bulk viscosity, respectively; and
$\boldsymbol{g} =-g\boldsymbol{e}_z$
is the gravitational acceleration. In terms of the above expansions, we will assume that the gravitational term does not appear until second order in
$M_a$
, which amounts to an assumption that
$\mu v_1/\ell ^2 \gg \rho g$
(where
$v_1$
is a representative value of
$|\boldsymbol{v}_1|$
and
$\ell$
is a representative droplet size), ensuring that gravity is negligible in both leading-order and first-order problems. This condition was checked a posteriori and found to hold. With these assumptions, we substitute the expansions from (3.1) into (3.2) and collect terms of like orders in
$M_a$
; §§ 3.1 and 3.2 discuss the first- and second-order terms, respectively.
3.1. First-order solution
The first order is associated with the propagation of acoustic waves in the fluid and hence we use the Helmholtz decomposition to separate the flow field into irrotational and divergence–free flow, i.e.
$\boldsymbol{v}_1=\boldsymbol{\nabla }\phi +\boldsymbol{\nabla }\times \varPsi$
, where
$\phi$
is the velocity potential of the irrotational flow component, and
$\varPsi$
represents the divergence–free flow component, so that the flow vorticity is equal to
$-{\nabla} ^2 \varPsi$
. We ignore contributions from
$\varPsi$
, assuming
$\boldsymbol{v}_1\approx \boldsymbol{\nabla }\phi$
. (The vorticity is known to be concentrated mainly in a boundary layer flow near the solid surface of thickness
$\sqrt {\mu /(\rho _0\omega )}$
, which amounts to just a few hundred nanometres and introduces only small contributions to the bulk flow far from the solid compared with the acoustic wave described by the potential
$\phi$
; see Vanneste & Bühler (Reference Vanneste and Bühler2011) for a general approach that further accounts for the boundary layer flow.) Here and in the following, it is understood that the real part of complex equations is taken when presenting results for real quantities such as
$\phi$
,
$\boldsymbol{v}_1$
or
$\boldsymbol{v}_2$
.
For simplicity of notation we then set
$\boldsymbol{v}_1 = \boldsymbol{\nabla }\phi$
in (3.1), leading (at order
$M_a$
on substitution in (3.2) and (3.3)) to the damped wave equation (Vanneste & Bühler Reference Vanneste and Bühler2011)
\begin{equation} \phi _{tt}=c^2{\nabla} ^2\phi +(\nu +\nu ' ){\nabla} ^2\phi _t, \end{equation}
where
$c^2=\left . {\rm d}p/{\rm d}\rho \right |_0$
,
$\nu '=\nu _b + \nu /3$
, with
$\nu \equiv \mu /\rho _0$
and
$\nu _b\equiv \mu _b/\rho _0$
. Assuming that the potential function,
$\phi$
, is harmonic in time with angular frequency
$\omega$
, and hence given by
$\phi =\hat {\phi }(\boldsymbol{x})\exp {(i\omega t)}$
, (3.4) may be reduced to
\begin{equation} {\nabla} ^2\hat {\phi }=- \frac {\omega ^2}{c^2+i(\nu +\nu ')\omega }\hat {\phi } = -\kappa ^2 \hat {\phi }, \end{equation}
where
\begin{equation} \kappa =\frac {k}{\sqrt {1+ i(\nu +\nu ') k/c}}\,, \end{equation}
and
$k= \omega /c=2\pi /\lambda$
is the wavenumber of the sound wave in the liquid, with wavelength
$\lambda$
. The real and imaginary parts of
$\kappa$
are given by
\begin{equation} \kappa = k_{{r}} - i \,k_{{i}}= k \sqrt {\cos \beta } \,e^{-i \beta /2} , \end{equation}
assuming
$\beta \in (0,( {\pi }/{2}))$
where
\begin{equation} \tan \beta = ( \nu +\nu ' ) \frac {k}{c}. \end{equation}
The solution of (3.5) requires a boundary condition that captures the corresponding leading-order behaviour for the normal component of a Rayleigh-type SAW at the solid surface. This condition can be derived from the solution given by Royer & Dieulesaint (Reference Royer and Dieulesaint1996) for the Rayleigh wave in an isotropic solid as
\begin{equation} \left . \frac {\partial \phi }{\partial z}\right |_{z=0}= A\omega \exp {(i\omega t)}\exp {(-i\kappa _{{s}} x)}, \end{equation}
where subscript
${s}$
stands for solid, and
$A$
is the displacement amplitude in the liquid. Unlike the solution in Royer & Dieulesaint (Reference Royer and Dieulesaint1996), the wavenumber
$\kappa _{{s}}$
has an imaginary part, i.e.
\begin{equation} \kappa _{{s}}=k_{{s},{r}}-i k_{{s},{i}}, \end{equation}
due to the presence of the liquid on top of the solid that converts the simple SAW into a LSAW. The consequent attenuation of the LSAW is accounted for by the imaginary part of
$\kappa _{{s}}$
, which is given by (Arzt, Salzmann & Dransfeld Reference Arzt, Salzmann and Dransfeld1967)
\begin{equation} k_{{s},{i}} = \frac {\rho _0}{\rho _{{s}}}\frac {c }{c_{{s}}^2}f, \end{equation}
where
$\rho _{{s}}$
is the solid density,
$f=\omega /2\pi$
is the SAW frequency and
$c_{{s}}$
(
$=\omega /k_{{s},{r}}$
) is the SAW phase velocity in the substrate. Note that the boundary condition, (3.9), is derived using the expansion of the nonlinear Rayleigh wave for small strains in the solid, so that the actual deformations at the solid surface as well as the inherently nonlinear relation between strain and stress in the solid, contribute at higher orders of
$M_a$
, not considered here.
We comment here also on the relation between the experimentally measured amplitude,
$A_{{n}}$
, and the theoretical one,
$A$
, used in (3.9). Although we expect
$A\propto A_{{n}}$
, it is not clear that these amplitudes are identical, as discussed in some detail by Royer & Dieulesaint (Reference Royer and Dieulesaint1996). A detailed analysis to determine the exact relation between
$A$
and
$A_{{n}}$
is beyond the scope of the present work, so we assume that
$A = A_{{n}}$
, a point to which we return later when discussing the comparison between experimental and theoretical results.
Values of the physical parameters related to the experiments and the derived lengths for both the solid substrate and the liquid (PDMS).
The experimental parameter values relevant to the silicone oil–lithium niobate system that we consider here are summarised in table 1. Note that, for the considered experimental parameters, we have
$\tan \beta \approx \beta \ll 1$
or
$\nu +\nu ' \ll c/k$
, so that to leading order in
$\beta$
(see (3.7), (3.8))
\begin{equation} k_{{r}} \approx k, \quad k_{{i}} \approx k \frac {\beta }{2} \approx \frac {k^2 \left ( \nu +\nu ' \right )}{2 c}. \end{equation}
This approximate value of
$k_{{i}}$
is in agreement with that defined in (5) of Brunet et al. (Reference Brunet, Baudoin, Matar and Zoueshtiagh2010). Note that
$\beta \ll 1$
for a range of fluids; e.g. for water (
$\nu =1$
cSt) and oil (
$\nu =50$
cSt) we obtain
$\beta =1.67\times 10^{-4}$
and
$\beta =8.38\times 10^{-3}$
, respectively, assuming
$\nu _b=\nu$
in both cases. Therefore, we can consider the ratio
$k_{{i}}/k_{{r}}$
as a small parameter, as discussed further in Appendix B.
The problem prescribed by (3.5) and (3.9), alongside a requirement that the acoustic field velocity decays far from the solid surface as
$z\rightarrow \infty$
, is satisfied by the potential function
\begin{equation} \phi =\frac {- A\omega }{\sqrt {\kappa _{{s}}^2-\kappa ^2}}\exp ({i\omega t})\exp {\left (-i\kappa _{{s}} x\right )} \exp {\left (-z\sqrt {\kappa _{{s}}^2-\kappa ^2} \right )}, \end{equation}
yielding the velocity components
\begin{align} v_{1,x}&=\frac {\partial \phi }{\partial x}=\frac {i\kappa _{{s}}}{\sqrt {\kappa _{{s}}^2-\kappa ^2}} A\omega \exp ({i\omega t}) \exp {\left (-i\kappa _{{s}} x\right )}\exp {\left (-z\sqrt {\kappa _{{s}}^2-\kappa ^2} \right )}, \nonumber \\ v_{1,z}&=\frac {\partial \phi }{\partial z}=A\omega \exp ({i\omega t})\exp {\left (-i\kappa _{{s}} x\right )} \exp {\left (-z\sqrt {\kappa _{{s}}^2-\kappa ^2} \right )}. \end{align}
For reference, the (physically relevant) real parts of these velocity components are written out explicitly in Appendix A.
We note that the analysis presented so far follows parts of the work by Campbell & Jones (Reference Campbell and Jones1970) and Vanneste & Bühler (Reference Vanneste and Bühler2011), albeit instead of calculating the SAW in the solid from conservation equations, we take the simpler approach of representing the component of the LSAW at the solid surface as a harmonic function that decays exponentially along its path (Arzt et al. Reference Arzt, Salzmann and Dransfeld1967) and hence are able to obtain the analytical representation of the acoustic contribution to the velocity field in the fluid. The formulation that ignores viscous dissipation is discussed in Appendix B.
3.2. Second-order solution
Having obtained
$\boldsymbol{v}_1$
, we next seek the quasi–steady flow field at times long compared with the acoustic period (Stokes Reference Stokes1847; Rayleigh Reference Rayleigh1884; Schlichting Reference Schlichting1932). Since we excite the fluid using a single-frequency SAW, it is sufficient to time average the system of equations over the period of the SAW (Vanneste & Bühler Reference Vanneste and Bühler2011) using the operator
$\langle \boldsymbol{\cdot }\rangle \equiv \omega /(2\pi ) \int _{t=0}^{2\pi /\omega }\boldsymbol{\cdot }\, {\rm d}t$
. The leading-order velocity contribution that supports a steady flow component (and hence does not vanish) is
$\langle \boldsymbol{v}_2\rangle$
, appearing at order
$M_a^2$
in the expansions specified by (3.1). The corresponding conservation of momentum and mass equations are then given by
\begin{align} -\boldsymbol{\nabla }\left \langle p_2\right \rangle + {\boldsymbol{F}}_{\!{s}} &-\rho _0 g\boldsymbol{e}_z + \mu {\nabla} ^2 \left \langle \boldsymbol{v}_2\right \rangle + \left ( \mu _b + \frac {\mu }{3} \right ) \boldsymbol{\nabla }\left ( \boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle \right )=0, \end{align}
\begin{align} &\quad \rho _0\boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle +\boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \rho _1\boldsymbol{v}_1\right \rangle =0, \end{align}
where
\begin{equation} \boldsymbol{F}_{\!{s}}= -\langle \rho _0 \left ( \boldsymbol{v}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{v}_1 + \boldsymbol{v}_1 \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \rho _0 \boldsymbol{v}_1 \right ) \rangle , \end{equation}
represents the acoustic force on the induced steady streaming flow. Moreover, considering the component of the continuity equation proportional to
$M_a$
,
$\partial \rho _1/\partial t+\rho _0\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_1=0$
and the solution for
$\boldsymbol{v}_1$
specified by (3.14), we see that the product
$\rho _1\boldsymbol{v}_1$
leads to an odd function over the SAW period, which results in
$\langle \rho _1\boldsymbol{v}_1\rangle =0$
. Hence, (3.15) and (3.16) may be simplified to
\begin{align} -\boldsymbol{\nabla }\left \langle p_2\right \rangle + {\boldsymbol{F}}_{\!{s}} &-\rho _0 g\boldsymbol{e}_z + \mu {\nabla} ^2 \left \langle \boldsymbol{v}_2\right \rangle =0, \end{align}
\begin{align} &\boldsymbol{\nabla }\boldsymbol{\cdot }\left \langle \boldsymbol{v}_2\right \rangle =0, \end{align}
with
${\boldsymbol{F}}_{\!{s}}$
as specified in (3.17). In Appendix A we use the real part of the first-order velocity field,
$\boldsymbol{v}_1$
, to reveal the form of the SAW forcing as given in (3.17) above, and in Appendix B we discuss a simplified version as well as the inviscid limit. The main points we emphasise are that the acoustic streaming force, given by (3.17), is non-conservative in contrast to the result in Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989), and that its prefactor is significantly dependent on the presence of viscous dissipation. The discussion is made explicit in Appendix A where the real component is shown along with the exact form of the coefficients therein.
3.3. Long-wave approximation
We carry out a standard long–wave expansion of (3.18) and (3.19), focusing on the simplified two-dimensional geometry that assumes translational invariance in the transverse
$y$
-direction (
$\partial /\partial y=0$
). The discussion of three–dimensional effects will be presented elsewhere. Here, we obtain the pressure,
$p \equiv \langle p_2\rangle$
, by substituting the
$z$
-component of
${\boldsymbol{F}}_{\!{s}}$
, given in (A7) and (A8), in the
$z$
-component of (3.18). We then integrate the resulting equation over the film thickness with respect to
$z$
, using the Laplace pressure boundary condition
$p(x,h(x,t))=-\gamma h^{\prime \prime }$
(simplifying the curvature of the free surface in the spirit of the long–wave approximation, and using
$\prime$
to denote
$\partial /\partial x$
), to find
\begin{equation} p(x,z) = -\gamma h^{\prime \prime } + \rho _0 g (h-z) +\frac {C_z P_0}{2K_z}\left [\psi (x,h)-\psi (x,z) \right ] , \end{equation}
where parameters
$K_z, C_z$
are specified in (A6), (A10) respectively, and we defined the auxiliary function
$\psi$
and parameter
$P_0$
as
\begin{equation} \psi (x,z) = e^{-2(k_{{s},{i}}x+K_zz)}, \qquad P_0=\rho _0 A^2 \omega ^2. \end{equation}
The last term on the right-hand side of (3.20) is the contribution of the bulk acoustic forcing to the local pressure in the film. To obtain the velocity component in the
$x$
-direction at this order,
$v_{2,x}$
, we substitute both the
$x$
-derivative of the pressure in (3.20), and
$F_{{s},x}$
as given in (A7), into the
$x$
-component of (3.18). We then integrate twice with respect to
$z$
, applying a zero shear-stress boundary condition at the free surface
$(\partial v_{2,x}/\partial z=0$
at
$z=h$
) and a no-slip condition at the substrate
$(v_{2,x}=0$
at
$z=0)$
, which yields
\begin{align} v_{2,x} &= -\frac {1}{2\mu } z(z-2h)\left [\gamma h^{\prime \prime \prime }-\rho _0 g h^{\prime } + C_z P_0\psi (x,h)\left (\frac {k_{{s},{i}}}{K_z}+h^{\prime }\right ) \right ] \nonumber \\ &\quad - \mathcal{C} \frac {P_0}{4\mu K_z^3}e^{-2k_{{s},{i}} x}\left (e^{-2K_zz}+2K_ze^{-2K_zh}z-1 \right ) , \end{align}
where
\begin{equation} \mathcal{C} = K_z C_x -k_{{s, i}} C_z , \end{equation}
with parameter
$C_x$
defined in (A9). Note that the last term in (3.22) leads to a velocity profile
$v_{2,x}$
that is modified from the usual parabolic shape in the neighbourhood of
$z=0$
. Interestingly, there is a critical value of
$\mathcal{C}$
, namely
$\mathcal{C}_{{crit}}\lt 0$
, such that for
$\mathcal{C}\lt \mathcal{C}_{{crit}}$
, the velocity profile includes negative values of
$v_{2,x}$
in the region near the substrate. We have not, however, observed such flow inversion with the present choice of parameters relevant to our experiments.
From the velocity profile above, (3.22), we obtain the film-averaged velocity and flux as
\begin{equation} u = \frac {1}{h}\int \limits _0^h v_{2,x} \ {\rm d}z, \quad Q=uh. \end{equation}
Conservation of mass for the fluid then requires that
\begin{equation} \frac {\partial h}{\partial t} + \frac {\partial Q}{\partial x}=0. \end{equation}
Upon using the velocity profile as obtained in (3.22) above, we find
\begin{equation} Q = -\frac {h^3}{3\mu }\frac {\partial {\mathcal P}}{\partial x} - \mathcal{C}\frac {P_0}{8\mu K_z^4}\psi (x,h)\left [ 2K_z^2h^2-1+e^{2K_zh}\left ( 1-2K_zh \right )\right ]\!, \end{equation}
where we have defined the effective pressure
\begin{equation} {\mathcal P} = -\gamma h^{\prime \prime } + \rho _0 gh+\frac {C_z P_0}{2K_z}\psi (x,h). \end{equation}
The first term on the right-hand side of
$\mathcal{P}$
is the capillary contribution to the quasi–steady pressure in the fluid film, wherein
$\gamma$
is the surface tension at the fluid–air interface, assumed constant; the second and third terms account for the gravitational body force; and the final term models the contribution from the acoustic force due to the SAW.
For our simulations, we non-dimensionalise the problem using an arbitrary length scale,
$\ell$
, a time scale
$t_c=3\mu \ell /\gamma$
such that leading-order terms balance in (3.25), and a pressure scale
$p_{{c}} =\gamma /\ell$
based on the capillary contribution in (3.27), so that we have
\begin{align} &\qquad (x,h)=\ell (\tilde x,\tilde h), \quad t =\frac {3\mu \ell }{\gamma } \tilde t, \quad {\mathcal P} = \frac {\gamma }{\ell } \tilde {{\mathcal P}}, \nonumber\\& (k_{{s},{i}},K_z,C_x,C_z)=\ell ^{-1} (\tilde {k}_{{s},{i}},\tilde {K}_z,\tilde {C}_x,\tilde {C}_z), {\mathcal C}=\ell ^{-2} \tilde {\mathcal C}. \end{align}
Therefore, in dimensionless form, we have
\begin{equation} \frac {\partial \tilde {h}}{\partial \tilde {t}}+\frac {\partial }{\partial \tilde {x}}\bigg [-\tilde {h}^3 \frac {\partial \tilde {{\mathcal P}}}{\partial \tilde {x}} - \tilde {\mathcal{C}} \frac {3\mathcal{S}}{8 \tilde {K}_z^4}\psi \big(\tilde {x},\tilde {h}\big)\big (2 \tilde {K}_z^2 \tilde {h}^2-1+e^{2 \tilde {K}_z \tilde {h}}\big(1-2 \tilde {K}_z \tilde {h}\big) \big ) \bigg]=0, \end{equation}
where it is understood that in calculating
$\psi (\tilde {x},\tilde {h})$
we use (3.21) with
$k_{ {s,i}}$
,
$K_z$
replaced by their tilded equivalents; the dimensionless effective pressure is
\begin{equation} \tilde {\mathcal P} = -\tilde {h}^{\prime \prime } + {\textit{Bo}}\, \tilde {h} + \frac {\mathcal{S} \tilde {C}_z}{2 \tilde {K}_z} \psi ( \tilde {x},\tilde {h}), \end{equation}
and we have defined the non-dimensional parameters
\begin{equation} {\textit{Bo}}\,=\frac {\rho _0 g \ell ^2}{\gamma }=\frac {\ell ^2}{a^2}, \quad \mathcal{S} = \frac {P_0 \ell }{\gamma } = \frac {\rho _0 \ell A^2 \omega ^2}{\gamma }, \end{equation}
where
$a=\sqrt {\gamma /(\rho _0 g)}$
is the capillary length. In our simulations, motivated by a typical drop size in the experiments, we set
$\ell = 1$
mm. Table 2 provides a list of non-dimensional parameter values used in all simulations unless otherwise stated.
Dimensionless parameters used in the simulation of the experiment for a drop volume of
$\mathcal{V}_{{d}}=8\, \unicode{x03BC}$
l driven by a SAW with
$A$
in the range from table 1, which determines the range of
$\mathcal{S}$
.
We pause here to provide a justification for neglecting the effects of acoustic radiation pressure in the present model. The experiments modelled here involve a SAW propagating in the solid with attenuation length of
$1$
–
$2$
mm, while the drops considered are
$6$
–
$7$
mm in lateral length. When comparing the relative strengths of acoustic streaming and acoustic radiation pressure, an important distinction is that streaming results from the SAW in the solid, while acoustic radiation pressure results from the LSAW in the liquid that reflects off the drop’s free surface. The attenuation lengths of both effects are proportional to the relative phase velocities of the acoustic fields in the solid and fluid, respectively. Since the speed of sound in the liquid is
$2$
–
$3$
times smaller than the speed of sound in the solid, the radiation pressure attenuates twice as fast as the SAW in the solid. Hence, the radiation pressure affects the drop rear only along a distance of approximately
$0.5$
–
$1$
mm. Because the present focus is on describing the front dynamics of driven drops, we neglect this contribution. However, we note that in cases where the attenuation length of the SAW is comparable to the lateral length of the drop, acoustic radiation pressure is expected to have an effect on the drop geometry and dynamics comparable to that of acoustic streaming, and hence should be included (see, e.g. Marcos et al. Reference Marcos, Li, Fasano, Diez, Cummings, Manor and Kondic2025). To quantify the regime in which radiation pressure becomes important, some insight can be reached by considering the non-dimensional ratio
\begin{equation} R=\frac {1}{2r_{{d}}k_{{i}}}, \end{equation}
that is, the attenuation length of the SAW in the solid divided by the lateral length of the drop. When
$R = \mathcal O (1)$
, acoustic radiation pressure should be considered to properly model the drop front, but when
$R\ll 1$
(e.g.
$R=1/5$
as in our case), it can be neglected when modelling the drop front dynamics.
4. Theoretical results
In this section, we discuss the predictions of the model developed in § 3, and compare them with the experimental results. In § 4.1, we present results for a few chosen reference cases, and then we discuss the influence of model parameters, including the acoustic amplitude,
$A$
, the droplet volume
$\mathcal{V}_{{d}}$
and the droplet kinematic viscosity
$\nu$
, on the results. In § 4.2 we discuss (approximate) travelling-wave solutions.
4.1. Main features of the results
In our simulations, we solve (3.29) numerically using COMSOL
$^{\textrm {TM}}$
; see Appendix C for the description of the implementation. As the initial condition, we consider the two–dimensional parabolic drop of (dimensionless) cross-sectional area
\begin{equation} \tilde {{\mathcal A}}_{{d}} = 2\int _0^{\tilde {r}_{{d}}} \tilde {h}_{{d}} \left [ 1- \left (\frac {\tilde {x}}{\tilde {r}_{{d}}}\right )^2 \right ] \,{\rm d} \tilde {x} = \frac {4}{3} \tilde {h}_{{d}} \tilde {r}_{{d}}, \end{equation}
(recall that tilded quantities are dimensionless) which remains constant during the drop evolution. The drop height parameter
$\tilde {h}_{{d}}$
is obtained from the known experimental (three-dimensional) drop volume,
${\mathcal V}_{{d}} =\ell ^3 \tilde {\mathcal V}_{{d}}$
, as
$\tilde {h}_{{d}} = 2 \tilde {\mathcal V}_{{d}}/(\pi \tilde {r}_{{d}}^2)$
, using the formula for the volume
$\tilde {\mathcal V}_{{d}}$
of a three-dimensional (3-D) parabolic cap of base radius
$\tilde {r}_{{d}}$
and height
$\tilde {h}_{{d}}$
. The 2-D parabolic initial condition is centred at
$\tilde {x}=\tilde {x}_{{d}}$
, and is written as
\begin{equation} \tilde {h}(\tilde {x},0)= \left \{ \begin{array}{ll} \tilde {h}_{{p}} ,& \tilde {x} \lt \tilde {x}_{{r}}(0) ,\\ \big( \tilde {h}_{{d}} - \tilde {h}_{{p}} \big) \left [ 1- \left (\frac { \tilde {x}- \tilde {x}_{{d}}}{ \tilde {r}_{{d}}}\right )^2 \right ] + \tilde {h}_{{p}} ,& \tilde {x}_{{r}}(0)\leq \tilde {x} \leq \tilde {x}_{{f}}(0),\\ \tilde {h}_{{p}} ,& \tilde {x} \gt \tilde {x}_{{f}}(0), \end{array} \right . \end{equation}
where
$\tilde {x}_{{r}}(0)= \tilde {x}_{{d}} - \tilde {r}_{{d}}$
and
$ \tilde {x}_{{f}}(0)= \tilde {x}_{{d}} + \tilde {r}_{{d}}$
are the initial positions of the rear and front contact lines of the drop, respectively (see figure 5). Here,
$ \tilde {h}_{{p}}$
(
$ \ll \tilde {h}_{{d}}$
) is the thickness of a precursor film, introduced to avoid the stress singularity associated with a moving contact line. It is well known (see, e.g. Diez & Kondic Reference Diez and Kondic2002), that the step size
$\Delta \tilde {x}$
in the spatial discretisation should be similar to
$\tilde {h}_{{p}}$
to ensure numerical convergence, and therefore we use
$\Delta \tilde {x} = \tilde {h}_{{p}}$
in our simulations. The specific value assigned to
$\tilde {h}_{{p}}$
has only a minor influence on the results; e.g. if
$\tilde {h}_{{p}}$
is doubled or halved, the change of spreading speed or drop height is just 1 %–2 %. All reported results are obtained using
$\tilde {h}_{{p}} = 0.01$
.
In implementing our model for the geometry defined above, we note that in the experiments, the SAW experiences minimal attenuation in regions of the actuator where fluid is not present, with significant attenuation only under the bulk droplet. To model this, we modify the factor
$\psi$
defined in (3.21), setting it to a constant value (no attenuation) in regions where the film height is smaller than some critical thickness
$\tilde {h}^*$
, with the exponential attenuation only where
$\tilde {h}\gt \tilde {h}^*$
, so that
\begin{equation} \psi (\tilde {x},\tilde {h}) = \left \{ \begin{array}{ll} 1, & \tilde {x}\lt \tilde {x}_1^\ast (\tilde {t})\, , \\ e^{-2 \left [ \tilde {k}_{{s},{i}} (\tilde {x} - \tilde {x}_1^\ast (\tilde {t})) + \tilde {K}_z \left ( \tilde {h} - \tilde {h}^\ast \right ) \right ]}, & \tilde {x}_1^\ast (\tilde {t}) \leq \tilde {x} \leq \tilde {x}_2^\ast (\tilde {t})\, , \\ e^{-2 \tilde {k}_{{s},{i}} \left [ \tilde {x}_2^\ast (\tilde {t}) - \tilde {x}_1^\ast (\tilde {t}) \right ]}, & \tilde {x} \gt \tilde {x}_2^\ast (\tilde {t}) \, , \end{array} \right . \end{equation}
where
$\tilde {x}_1^\ast$
,
$\tilde {x}_2^\ast$
are defined by
$\tilde {h}(\tilde {x}_1^\ast ,\tilde {t})=\tilde {h}(\tilde {x}_2^\ast ,\tilde {t})=\tilde {h}^\ast$
, and the SAW is assumed to propagate from
$\tilde {x}\lt \tilde {x}_1^*(\tilde {t})$
in the direction of increasing
$\tilde {x}$
. In line with our discussion in § 3, in simulations we choose the cutoff film thickness
$\tilde {h}^*$
based on the thickness below which the bulk acoustic streaming becomes weak and mass transport is governed by other mechanisms (e.g. Rayleigh streaming and acoustic pressure Rezk et al. Reference Rezk, Manor, Friend and Yeo2012, Reference Rezk, Manor, Yeo and Friend2014),
$\tilde {h}^\ast \approx (\lambda _{{oil}}/4)/\ell =17 \,\unicode{x03BC} \textrm {m}/\ell$
.
We commence our analysis by considering the experimental case where an oil drop of volume
$8\, \unicode{x03BC}$
l is under the excitation of the
$20$
MHz SAW. Figure 6 shows the evolution of the thickness profile obtained by solving (3.29) for
$A=1.44, 2, 6, 10$
nm subject to the initial condition given by (4.2), and boundary conditions fixing the film thickness to
$h_{{p}}$
and requiring vanishing derivatives at the domain boundaries; the other parameter values are shown in table 2. For the low-amplitude cases (e.g.
$A\lt 2$
nm, figures 6
a and 6
b), we see that capillary and gravitational forces dominate, leading to a continual spreading of the drop. The acoustic force is able to slightly drive the drop along the direction of propagation of the SAW, as is apparent by the movement of the two contact lines. In contrast, in the high-amplitude cases (e.g.
$A\gt 2$
nm, figures 6
c and 6
d), the acoustic force plays a dominant role in governing the dynamics. We see that the evolution of the drop can be described by two regimes: an initial transient stage where the drop takes on a new shape, followed by a period of translation across the domain at near-constant speed, consistent with experimental results (see figure 3); as
$A$
increases, the duration of the transient stage decreases. The fluid region near the rear contact line, where the pressure gradient is largest, serves as a ‘snow plough’ that pushes the bulk drop along the direction of propagation of the SAW. The result is a near trave-ling-wave solution at long times, discussed further in § 4.2. The subtle, but continuous, decrease of the maximum height of the drop during this period is due to a thin trailing film (of approximate thickness
$h^\ast$
) that is left behind the bulk drop. We also note formation of a trailing foot forming spontaneously behind the main body of the drop for amplitudes
$A\gt 2$
nm. This feature travels together with the main body of the droplet, in contrast to the trailing film, which stays behind. We direct the interested reader to the supplementary materials where animations for each of these four cases are shown (movie 4, movie 5, movie 6, movie 7).
Evolution of thickness profile for different values of
$A$
. The black solid lines correspond to the initial condition, the dashed black lines to the profiles every time step
$\Delta t=1$
s and the solid blue line to
$t=6\,$
s. The red line stands for the effective dimensionless pressure
$\tilde {\mathcal{P}}$
at
$t=6$
s. Note that the dip in the pressure curve corresponds to the dramatic change of curvature of the fluid interface at the front contact line.
Focusing next on the spreading speed, figure 7 plots the relative speed of spreading for both experiments and simulations. The numerical results (solid lines) show a monotonic increase in spreading speed as
$A$
increases, consistent with the experimental results (symbols). Note the distinct regimes in which the drops evolve, evident from the apparent change of slope of the curves, especially prominent in large-amplitude cases. Figure 7(b) shows the numerically calculated drop width,
${w}=x_2^\ast -x_1^\ast$
, as a function of time, for different values of
$A$
. We once again emphasise the different behaviours exhibited by low- and high-amplitude solutions. For low-amplitude cases, the front contact line of the drop moves faster than the rear contact line, causing the drop to become thinner and wider. Conversely, for the high-amplitude cases,
$\rm {w}$
initially increases and then remains nearly constant, or even decreases slowly.
Figure 7(a) illustrates a quantitative difference between experiment and simulations: assuming that
$A = A_{{n}}$
(see also the relevant discussion in § 3), the simulated drops travel at a slower pace. While it is not yet clear why this difference arises, some plausible reasons, in addition to possible differences between
$A$
and
$A_{{n}}$
(Royer & Dieulesaint Reference Royer and Dieulesaint1996), include: (i) the potential relevance of inertial effects that are not included in the model at order
$M_a$
(the solution
$\boldsymbol{v}_1$
) (Zarembo Reference Zarembo1971; Orosco & Friend Reference Orosco and Friend2022; Dubrovski, Friend & Manor Reference Dubrovski, Friend and Manor2023); (ii) the omission of contributions from the boundary layer flow near the solid surface, i.e. the Rayleigh law of streaming; (iii) the capillary waves that appear at the free surface of the film (Lighthill Reference Lighthill1978; Morozov & Manor Reference Morozov and Manor2017); (iv) the use of a 2-D model to understand a 3-D system; and (v) the neglect of the ultrasonic wave reflection off the free surface of the liquid. We have verified that the numerical results are independent of the precise values assigned to
$h^*$
and
$h_{{p}}$
, suggesting that neither boundary layer thickness nor precursor thickness are relevant here. We have also carried out some preliminary simulations of a 3-D oil drop, which suggest that the 2-D model used in the present work is not responsible for the quantitative differences between experiments and simulations (we leave the discussion of 3-D effects for future work). Clearly, more work will be needed to uncover the sources of the quantitative differences between experiments and simulations.
(a) Comparison between numerical (solid lines, amplitude values
$A$
as labelled, colour coded) and experimental results (solid dots, for measured amplitude values
$A_{{n}}=1.60\,$
nm (blue),
$1.44\,$
nm (red),
$0.95\,$
nm ( olive)). (b) Numerically calculated drop width data plotted as a function of time for several
$A$
-values (colour matched with panel (a) where appropriate).
Long time (calculated at
$t=9$
s) values of the front speed,
${\textit{v}}_{{f}}$
(in mm s−1, blue filled circles), maximum thickness,
$h_{\textit{max}}$
(mm, red filled squares) and drop width,
$\rm w$
(cm, black filled circles), as functions of
$A$
. Note that
$\textit{v}_{{f}}$
is plotted using the left
$y$
-axis while
$\rm w$
and
$h_{\textit{max}}$
are plotted using the right
$y$
-axis. The points correspond to the raw data, and the solid lines connecting them guide the eye.
Evolution of thickness profile for
$A=2$
nm for drops of different viscosities,
$\nu =\nu '=50,100,500$
cSt, indicated in each panel. The black solid lines correspond to the initial condition, the dashed black lines to the profiles every time step
$\Delta t=2$
s and the solid blue line to
$t=10\,$
s. The red line denotes the effective dimensionless pressure
$\tilde {\mathcal{P}}$
at
$t=10$
s.
To summarise the numerical results presented so far, figure 8 plots the numerically predicted values of the front speed,
$\textit{v}_{{f}}$
, the maximum drop height,
$h_{\textit{max}}$
and the drop width,
$\rm {w}$
, all calculated at
$t=9$
s, as a function of SAW amplitude
$A$
. This figure illustrates the transition between the two types of dynamics observed: for low and moderate
$A$
-values (
$A\leq 2$
nm) we observe near-symmetric drop spreading due primarily to capillary forces, however, for large values (
$A\gt 2$
nm) the front speed increases approximately quadratically as
$A$
increases, while the maximum height/width of the drop increases/decreases approximately linearly with respect to the amplitude value.
Before concluding this section, we briefly discuss how the results depend on the oil viscosity. We have also carried out simulations using different drop volumes; while the results are consistent with the experimental ones, shown in figure 4(b), the influence of the change of fluid volume is, however, weaker than in the experiments (the results not shown for brevity). When varying the kinematic viscosity,
$\nu$
(assuming
$\nu '=\nu$
), without acoustic forcing (
$A=0$
), viscosity only (linearly) affects the time scale. With acoustics, however, the dynamics changes in a nonlinear fashion, see figure 9. The mechanism behind this nonlinear dependence is traced to the prefactor of the non-conservative term in (3.29), which decreases as
$\tilde {K}_z$
(defined in (A6)) grows and
$\tilde {\mathcal{C}}$
(defined in (3.23)) decreases in magnitude, yielding a relatively large attenuation length. This substantially weakens the forcing at the front, while nearly balancing gravitational and capillary forces at the rear, effectively stalling the drop motion for large viscosity values, in qualitative agreement with the experimental findings of figure 4(a).
4.2. Travelling-wave solution
In view of the results presented so far, we now briefly explore a possible trave-ling-wave solution of (3.29) which appears present for sufficiently large SAW amplitudes, such as the case shown in figures 6(c) and 6(d).
Within this framework, we assume that the whole drop translates with constant speed
$\tilde {U}$
, so that
\begin{equation} \tilde {x}_{{r}}(\tilde {t}) = \tilde {U} \tilde {t}, \quad \tilde {x}_{{f}}(\tilde {t}) = \tilde {\rm {w}} + \tilde {U}\tilde {t}, \end{equation}
where (constant)
$\tilde {\rm {w}}$
is the width of the moving drop. Here, we will consider that the rear (front) contact lines of the travelling drop are given by
$\tilde {x}_1^\ast$
(
$\tilde {x}_2^\ast$
) (see figure 5), since the SAW force that drives it is felt only where
$\tilde {h}\gt \tilde {h}^\ast$
. Then, we define
$\tilde {\xi } = \tilde {x} - \tilde {x}_1^\ast ( \tilde {t} )= \tilde {x} - \tilde {U} \tilde {t}$
and assume that
$\tilde {h}(\tilde {x},\tilde {t})=\tilde {h}(\tilde {\xi } )$
, so that (3.29) becomes (see also Buckingham, Shearer & Bertozzi Reference Buckingham, Shearer and Bertozzi2003; Perazzo & Gratton Reference Perazzo and Gratton2004; Giacomelli, Gnann & Otto Reference Giacomelli, Gnann and Otto2016 for similar studies without SAW forces)
\begin{align} &\tilde {U} \frac {{\rm d} \tilde {h}}{{\rm d} \tilde {\xi }} + \frac {\rm d}{{\rm d} \tilde {\xi }} \left [ \tilde {h}^3 \frac {d}{{\rm d} \tilde {\xi }}\left ( -\frac {{\rm d}^2 \tilde {h}}{{\rm d} {\tilde {\xi }}^2} + {\textit{Bo}} \, \tilde {h} + \frac {\mathcal{S} \tilde {C}_z}{2 \tilde {K}_z} \psi (\tilde {\xi } ,\tilde {h}) \right )\right . \nonumber\\& \quad \left . +\, \tilde {\mathcal{C}} \frac {3\mathcal{S}}{8 \tilde {K}_z^4}\psi ( \tilde {\xi },\tilde {h}) \bigg (2 \tilde {K}_z^2 \tilde {h}^2-1+e^{2 \tilde {K}_z \tilde {h}}(1-2 \tilde {K}_z \tilde {h}) \bigg )\right ] = 0, \end{align}
where
\begin{equation} \psi (\tilde {\xi } ,\tilde {h})= e^{-2 [ \tilde {k}_{{s},{i}} \tilde {\xi } + \tilde {K}_z (\tilde {h}- \tilde {h}^\ast )]}. \end{equation}
This equation can be integrated once to yield
\begin{align} &\tilde {U} \tilde {h} + \tilde {h}^3 \frac {\rm d}{{\rm d} \tilde {\xi }} \left [ -\frac {{\rm d}^2 \tilde {h}}{{\rm d} {\tilde {\xi }}^2} + {\textit{Bo}}\, \,\tilde {h} + \frac {\mathcal{S} \tilde {C}_z}{2\tilde {K}_z} \psi (\tilde {\xi }, \tilde {h}) \right ] \nonumber\\& \quad + \tilde {\mathcal{C}} \frac {3\mathcal{S}}{8 \tilde {K}_z^4} \psi (\tilde {\xi },\tilde {h}) \bigg (2 \tilde {K}_z^2 \tilde {h}^2-1+e^{2 \tilde {K}_z \tilde {h}}(1-2 \tilde {K}_z \tilde {h})\bigg )= \tilde {J}, \end{align}
where
$\tilde {J}$
represents the flux. The travelling-wave solution must be calculated for
$0 \leq \tilde {\xi } \leq \tilde {{w}}$
, along with the following boundary conditions at
$\tilde {\xi } =0$
:
\begin{equation} \tilde {h} (0)= \tilde {h}^\ast , \quad \tilde {h}^\prime (0)= \tilde {h}^{\prime \prime \prime }(0)=0, \end{equation}
where the prime denotes
$d/{\rm d} \tilde {\xi }$
. The values of
$\tilde {h}^{\prime \prime }(0)$
,
$\tilde {U}$
and
$\tilde {{w}}$
are determined by
\begin{equation} \tilde {h} (\tilde {{w}})= \tilde {h}^\ast , \quad \tilde {J}(0)= \tilde {J} (\tilde {{w}}), \quad \tilde {A}_{{d}}=\int _0^{\tilde { {w}}} \tilde {h}(\tilde {\xi })\,{\rm d}\tilde {\xi }. \end{equation}
Note that
$\tilde {J}$
corresponds to the flux within a thin film of thickness
$\tilde {h}^\ast$
that enters the drop at
$\tilde {{w}}$
and comes out of it at
$\tilde {\xi } =0$
(with respect to the reference frame fixed at the drop).
We developed an iterative scheme to perform the numerical integration of (4.7), under the conditions given by (4.9). We start by guessing values of
$\tilde {h}^{\prime \prime }(0)$
,
$\tilde {U}$
and
$\tilde {\xi }_{{f}}$
to perform the integration of (4.7) and modify them accordingly until these conditions are satisfied within a small relative error (typically,
$10^{-7}$
or even smaller). After obtaining convergence for a given
$\tilde {A}$
, we use these converged values of
$\tilde {h}^{\prime \prime }(0)$
,
$\tilde {U}$
and
$\tilde {\xi }_{{f}}$
as new guess values for
$\tilde {A} + \Delta \tilde {A}$
(numerical continuation;
$\Delta \tilde {A}$
has to be very small to reach convergence).
Figure 10(a) shows the dimensional drop profiles for the values of
$A$
considered so far. Clearly, drops become taller and narrower as
$A$
increases. Figure 10(b) shows the dependence of drop speed,
$\textit{v}_{{f}}$
, its width,
$\rm w$
, and the maximum thickness,
$\textit{h}_{\textit{max}}$
, on
$A$
. Note that this simplified model predicts that the speed depends approximately quadratically on the amplitude,
$A$
, while the maximum height of the drop depends approximately linearly on
$A$
, as was found in the full numerics (see figure 8
b).
This approximate formulation shows that the simple travelling-wave solution can capture, at least qualitatively, certain features of the experimental and numerical drops. Clearly, this description cannot account for the rear region of the drop (see figures 6 and 10
a) because the values of
$h'$
and
$h'''$
at the left contact line cannot be determined a priori, and both have been set to zero (for lack of better choice) in the travelling-wave calculations.
(a) Thickness profile of the travelling-wave solution for:
$A=4.5$
nm,
$6$
nm,
$8$
nm,
$10$
nm. (b) Front speed,
$\textit{v}_{{f}}$
(mm s−1, blue line), drop width,
$\rm w$
(mm, black line) and maximum height,
$h_{{max}}$
(mm, red line), of the travelling drop as a function of
$A$
; front speed is plotted on the left
$y$
-axis while drop width and maximum height are on the right
$y$
-axis.
5. Conclusions
In this work, we consider the problem of a spreading silicone oil drop driven by a SAW propagating in the supporting substrate. Our particular focus is on a theoretical framework that describes drops of millimetric thickness, for which the main driving mechanism is the Eckart streaming. Therefore, we focus on a different regime from that relevant for the much thinner drops considered by Rezk et al. (Reference Rezk, Manor, Friend and Yeo2012, Reference Rezk, Manor, Yeo and Friend2014), for which Schlichting and Rayleigh streaming effects are dominant. To the best of our knowledge, modelling the dynamics of millimetric drops under the action of SAW and capillary stresses has not been considered so far in the literature. We also present selected experimental results for the qualitative validation of our model.
While developing our theoretical description of drop spreading, it became apparent that an important effect, which requires careful modelling, involves the attenuation of the SAW under the evolving drop. The model we derived, assuming the thin film limit and implementing the long-wave approximation, shows that the oil drop is essentially pushed by the SAW forcing from behind. For the highly viscous oils considered in our experiments, the effective SAW forcing is due to both viscous dissipation, and to spatial variation of the leaked acoustic wave in the fluid due to SAW attenuation in the solid. A consequence of the driving mechanism is that the dynamics itself is substantially different from the familiar case of drops driven by gravity, where the volume force is a constant and not rapidly decreasing along the drop, as occurs with the SAW force. In particular, we find that the drops tend to keep almost the same shape as time and spreading progress, and also that they reach essentially constant spreading speeds at long times. This contrasts with the gravity-driven case, where the spreading speed decreases due to drop thinning.
While not all details of the experimental results are fully understood, in general, we find reasonable agreement between the experimental and theoretical results, further supporting the basic premises underlying the developed acousto-fluidic model. In particular, our model provides good qualitative predictions for drop shapes, spreading speeds and the dependence of the results on the SAW displacement amplitude at the solid surface (the SAW intensity). Both theoretical and experimental results show asymmetric drop shapes with a thin trailing film left behind, and both appear to reach a constant spreading speed asymptotically in time. The latter finding motivated us to formulate a simplified travelling-wave model. The resulting analytical solution predicts an approximately linear increase in the drop thickness and a quadratic increase of the spreading speed with the acoustic intensity (SAW normal displacement amplitude at the solid surface),
$A$
. Such predicted trends are consistent with those of the full theoretical model and the experiments.
While many features of the experimental results have been rationalised using our theoretical models and simulations, much remains to be done. Further investigation is required to shed light on the features of dynamic drop profiles under SAW excitation. The transition from a 3-D to a quasi-2-D dynamics is also intriguing: under strong SAW excitation, a drop spreads primarily along the path of the SAW and may be described to leading order using the 2-D model given here. However, at weak SAW excitation, where the drop also spreads radially due to the relaxation of capillary stress at the free drop surface, the drop dynamics is inherently of the 3-D type. Hence, 3-D aspects of the drop dynamics, in particular in cases where the contributions to flow from capillary and SAW stresses are comparable, remain to be simulated and understood.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10802.
Acknowledgements
This work was supported by the donors of ACS Petroleum Research Fund under PRF# 62062-ND9 and by BSF grant No. 2020174. J.A.D. acknowledges support from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina) with Grant PIP 02114-CO/2021 and Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, Argentina) with Grant PICT 02119/2020.
Declaration of interests
The authors report no conflict of interest.
Appendix A. First-order velocity components and SAW force
As mentioned in the main text, although we deal with imaginary equations in the model derivation, it is understood that the real part of these complex equations is taken as the result for physically relevant real quantities. This includes the first-order velocity field,
$\boldsymbol{v}_1$
, given by (3.14), and the acoustic force on the induced steady streaming flow,
$\boldsymbol{F}_{\!{s}}$
, given by (3.17), which we now discuss in detail.
We first examine
$\boldsymbol{v}_1$
whose real parts in two dimensions can be written as
\begin{align} v_{1,x} &= V_xA\omega e^{-k_{{s},{i}} x-K_zz}, \nonumber \\ v_{1,z} &= V_zA\omega e^{-k_{{s},{i}} x-K_zz}, \end{align}
where the factors
$V_x$
and
$V_z$
are
\begin{align} V_x &= \frac {1}{\Delta } \left [ k_{{s},{i}} \cos \left (\frac {\theta }{2}+k_{{s},{r}} x+ z \Delta \sin \left (\frac {\theta }{2}\right )- \omega t \right )\right .\nonumber\\&\quad\left . +\,k_{{s},{r}} \sin \left (\frac {\theta }{2}+k_{{s},{r}} x+ z\Delta \sin \left (\frac {\theta }{2} \right )-\omega t \right ) \right ]\!, \end{align}
\begin{align} V_z &= \cos \left (k_{{s},{r}} x + z \Delta \sin \left (\frac {\theta }{2}\right )-\omega t \right )\!, \end{align}
with
\begin{align} \Delta &= \sqrt {\left | \kappa _{{s}}^2 - \kappa ^2 \right |}=\sqrt [4]{\left (k_{{s},{i}}^2-k_{{s},{r}}^2+k_{{r}}^2-k_{{i}}^2\right )^2+(2 k_{{s},{i}} k_{{s},{r}}-2 k_{{i}} k_{{r}})^2}, \end{align}
\begin{align} \theta &= \arg \left ( \kappa _{{s}}^2 - \kappa ^2 \right )\!, \end{align}
\begin{align} K_z &= \Delta \cos \left ( \theta / 2 \right ). \end{align}
Turning next to
$\boldsymbol{F}_{\!{s}}$
, in order to find the real part we substitute (A1) with the subsequently defined parameters into (3.17), and find explicit expressions for the acoustic force components in the form
\begin{align} F_{{s},x} &= C_x P_0 \, e^{-2 (k_{{s},{i}} x + K_z z)}, \end{align}
\begin{align} F_{{s},z} &= C_z P_0 \, e^{-2 (k_{{s},{i}} x + K_z z)}, \end{align}
where
$P_0$
is defined in (3.21), and
\begin{align} C_x &= \frac {k_{{s},{i}} k_{{s},{r}}^2+\Delta \cos \left (\dfrac {\theta }{2}\right ) k_{{s},{i}} K_z+k_{{s},{i}}^3+\Delta \sin \left (\dfrac {\theta }{2}\right ) k_{{s},{r}} K_z}{\Delta ^2}, \end{align}
\begin{align} C_z &= \frac {\cos \left (\dfrac {\theta }{2}\right ) k_{{s},{i}}^2+\sin \left (\dfrac {\theta }{2}\right ) k_{{s},{i}} k_{{s},{r}}+\Delta \, K_z}{\varDelta }. \end{align}
We point out that
$\boldsymbol{F}_{\!{s}}$
is a non-conservative force; it cannot be written as the gradient of a potential, except for a special choice of parameter values that satisfy
$k_{{s},{i}} C_z = K_z C_x$
(see (A7) and (A8)).
Appendix B. Approximate solution of the SAW force
In this section, we discuss an approximate solution for the SAW force in certain contexts. Section B.1 discusses simplifications to the model considered in Appendix A, using the appropriate parameters for Polydimethylsiloxane (PDMS), as considered in the experiments and in the main body of the paper. To illustrate the differences that occur in the case of less viscous fluid (water), in § B.2 we discuss appropriate simplifications for such cases. Finally, in § B.3 we discuss the model that results if viscous losses are neglected, and the relation of the obtained model to that proposed decades ago by Shiokawa, Matsui & Ueda (Reference Shiokawa, Matsui and Ueda1990).
B.1. Oil (PDMS)
In this section, we take a closer look at the acoustic force
$\boldsymbol{F}_{\!{s}}$
, given by (3.17), on the induced steady streaming flow. We first note that, for our parameter values, we can define the small quantities
\begin{equation} \varepsilon _{{s}} = \frac {k_{{s},{i}}}{k_{{s},{r}}} \quad \text{and}\quad \varepsilon =\frac {k_{{i}}}{k_{{r}}}\approx \frac {\beta }{2}. \end{equation}
These two parameters
$\varepsilon _{{s}}$
and
$\varepsilon$
are always small, and of the same order of magnitude, for the PDMS used here (see (3.12) and the surrounding discussion). To approximate the SAW force, we can expand the auxiliary variables and SAW force coefficients, given in (A4), (A5), (A6), (A9) and (A10), in terms of
$\varepsilon _{{s}}$
and
$\varepsilon$
. Retaining only linear terms in the expansions (neglecting any terms quadratic in
$\varepsilon _{{s}}, \varepsilon$
that are of higher order than those retained) yields the following expressions:
\begin{eqnarray} \Delta &\approx & k_{{s},{r}} \alpha _1, \nonumber \\ \theta &\approx & \pi - \frac {2}{\alpha _1^2}\big(1+\alpha _1^2\big)\varepsilon + \frac {2}{\alpha _1^2}\varepsilon _{{s}}, \nonumber \\ K_z &\approx & \frac {k_{{s},{r}}}{\alpha _1}\big(1+\alpha _1^2\big)\varepsilon -\frac {k_{{s},{r}}}{\alpha _1}\varepsilon _{{s}}, \\ \nonumber C_x &\approx & \frac {k_{{s},{r}}\left(1+\alpha _1^2\right)}{\alpha _1^2}\varepsilon , \\ \nonumber C_z &\approx & \alpha _1 C_x, \end{eqnarray}
where we have simplified the expressions by defining
\begin{equation} \alpha _1 = \sqrt {\frac {k^2}{k_{{s},{r}}^2}-1} = \sqrt {\frac {c_{{s}}^2}{c^2}-1} , \end{equation}
since
$k_{{s},{r}} =\omega /c_{{s}}$
. This coefficient yields the Rayleigh angle
$\theta _R$
as:
$\tan \theta _R = 1/\alpha _1$
.
From the outlined approximation, we obtain that the resultant force is given by
\begin{align} F_{{s},x} &\approx \frac {\rho _0 A^2 \omega ^2 (1+\alpha _1)^2 k_{{s},{r}}}{\alpha _1^2}\varepsilon \ e^{-2(k_{{s},{i}}x+K_z z)}, \end{align}
\begin{align} F_{{s},z} &\approx \alpha _1 F_{{s},x} . \end{align}
The approximations made in (B2) lead to an attenuation factor
$K_z$
that is approximately
$12$
% smaller and force coefficients
$C_x$
and
$C_z$
that are approximately
$2$
% bigger and
$6$
% smaller, respectively, with respect to the quantities in the exact expressions given in (A6), (A9), and (A10). Note that the approximate force components obtained in (B4), (B5) (also) do not yield a conservative force, since
$K_z \neq \alpha _1 k_{{s},{i}}$
.
B.2. Water
Next, we consider modelling a lower viscosity fluid such as water, characterised by a similar density and phase velocity to oil (PDMS) so that the imaginary component of the wavenumber in the solid is similar (see (3.11)). However, the large decrease in viscosity in the case of water leads to a much reduced value of the imaginary wavenumber in the liquid,
$k_{{i}}$
(see (3.7)). Thus for water, our assumption that
$\varepsilon$
and
$\varepsilon _{{s}}$
are of the same order of magnitude fails; in fact, here we have
$\varepsilon _{{s}}^2 \sim \varepsilon$
so that we need to also retain terms that are quadratic in
$\varepsilon _{{s}}$
; the next largest terms in the expansion are those proportional to
$\varepsilon _{{s}}^3$
which we also include for the purpose of comparison with the approximation in the case of neglecting viscous dissipation which is discussed later in § B.3. Note that we neglect any cross-terms that may appear in the expansions as they are higher order than the terms retained.
Performing an asymptotic expansion of the auxiliary variables and SAW force coefficients appropriate for the water leads to the following approximate expressions:
\begin{align} \Delta &\approx k_{{s},{r}}\alpha _1 + \varepsilon _{{s}}^2 \bigg (\frac {k_{{s},{r}}}{\alpha _1^3}+\frac {k_{{s},{r}}}{2\alpha _1} \bigg ), \nonumber \\ \theta &\approx - \pi + \frac {2}{\alpha _1^2}\varepsilon _{{s}} - \frac {2\left(1+\alpha _1^2\right)}{\alpha _1^2}\varepsilon , \nonumber \\ K_z &\approx \frac {k_{{s},{r}}}{\alpha _1}\varepsilon _{{s}}-\frac {k_{{s},{r}}\left(1+\alpha _1^2\right)}{\alpha _1}\varepsilon , \\ \nonumber C_x &\approx \frac {k_{{s},{r}}\left(1+\alpha _1^2\right)}{\alpha _1^2}\bigg (\varepsilon +\frac {\varepsilon _{{s}}^3}{\alpha _1^2} \bigg ) , \\ \nonumber C_z &\approx \frac {k_{{s},{r}}\left(1+\alpha _1^2\right)} {\alpha _1}\bigg (-\varepsilon +\frac {\varepsilon _{{s}}^3}{\alpha _1^4} \bigg ). \end{align}
We pause to discuss the difference between the approximated complex arguments in the case of oil and water. Note that the two expressions for
$\theta$
((B2) and (B6)) have the same first-order correction terms with respect to the small parameters
$\varepsilon$
and
$\varepsilon _{{s}}$
; however, the sign of the sum of these two terms is opposite for the case of PDMS (negative) and water (positive) – which is due to the much reduced value of the imaginary wavenumber in the liquid,
$k_{{i}}$
, in the case of water as compared with PDMS (see (3.6) and (3.8)). The result of this sign change has the effect of moving the complex argument from quadrant II (in the case of PDMS) to quadrant III (in the case of water); as a result, the approximation for the sine function in the SAW force coefficients changes sign and thus the overall value of the SAW force coefficients is much reduced (see (A9) and (A10)).
From the approximation in (B6), we can see that the resultant force for water, at leading order, is given by
\begin{align} F_{{s},x} &\approx \frac {\rho _0 A^2 \omega ^2 (1+\alpha _1)^2 k_{{s},{r}}}{\alpha _1^2}\varepsilon \ e^{-2(k_{{s},{i}}x+K_z z)}, \end{align}
\begin{align} F_{{s},z} &\approx -\alpha _1 F_{{s},x} . \end{align}
Using appropriate parameter values for water, we find that the SAW force for water, (B7) and (B8), has a coefficient that is approximately two orders of magnitude smaller than that for PDMS, see (B4) and (B5) (due to the fact that
$\varepsilon$
is much smaller for water than oil). Additionally, the
$z$
component of the SAW force has a negative prefactor,
$C_z$
.
Before closing this discussion, we note that if one wants to accurately study the dynamics of water drops in the present context, one must also include a disjoining pressure term in (3.27) to account for the partial wettability of the fluid.
B.3. Neglecting viscous dissipation
We now discuss briefly a modified model obtained if viscous dissipation is neglected by setting
$k_{{i}}=0$
, which is akin to the approach presented by Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989, Reference Shiokawa, Matsui and Ueda1990). To do this, the simplest approach is to set
$k_{{i}}=0$
in the exact formulae ((A4), (A5), (A6), (A9), and (A10)) and then re-expand the quantities in terms of the remaining small parameter
$\varepsilon$
. Doing this, one finds that the auxiliary variables and the SAW force coefficients are now given by
\begin{eqnarray} \Delta &\approx & k_{{s},{r}} \alpha _1, \nonumber \\ \theta &\approx & -\pi + \frac {2}{\alpha _1^2}\varepsilon _{{s}}, \nonumber \\ K_z &\approx & \frac {k_{{s},{r}}}{\alpha _1}\varepsilon _{{s}} , \\ \nonumber C_x &\approx & \frac {k_{{s},{r}}\left(1+\alpha _1^2\right)}{\alpha _1^4}\varepsilon _{{s}}^3,\\ \nonumber C_z &\approx & \frac {1}{\alpha _1} C_x. \end{eqnarray}
We note that this approximation is exactly what is found by simply setting
$\varepsilon =0$
in the approximation for the case of a low viscosity fluid such as water (see (B6)). Hence, if viscous dissipation is neglected, the SAW force is given as
\begin{align} F _{{s},x} &\approx \frac {\rho A^2 \omega ^2 \left(1+\alpha _1^2\right)}{\alpha _1^4} k_{{s},{r}} \varepsilon _{{s}}^3 e^{-2k_{{s},{i}}(x+\frac {1}{\alpha _1} z)}, \end{align}
\begin{align} F_{{s},z} &\approx \frac {1}{\alpha _1} F_{{s},x} , \end{align}
which can be rewritten in terms of a potential
$\phi _{{s}}$
, so that
$\boldsymbol{F}_{\!{s}}=-\boldsymbol{\nabla }\phi _{{s}}$
is a conservative volume force with
\begin{equation} \phi _{{s}} = \frac {\rho A^2 \omega ^2 \left(1+\alpha _1^2\right)}{2\alpha _1^4 }\varepsilon _{{s}}^2 e^{-2k_{{s},{i}}(x+\frac {1}{\alpha _1}z)}\,. \end{equation}
We note that without viscous dissipation, the force components are five orders of magnitude smaller than when the viscosity is considered, compare for example (A9), (A10) to (B9) with the given parameter values. As a result, the force in (B10) and (B11) has no effect on drop dynamics on the time scales investigated. Such a finding inspires comparison with the results presented in Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989, Reference Shiokawa, Matsui and Ueda1990), and given in terms of our variables as
\begin{eqnarray} \phi _{\textrm {shio}} = \frac {1+\alpha _1^2}{2}\rho A^2 \omega ^2 e^{-2k_{{s},{i}}(x+\alpha _1 z)}. \end{eqnarray}
Note that, for the considered parameter values, this expression leads to a force that is approximately three orders of magnitude larger compared with what is found from our model, i.e. (B12). Further, our model predicts a different attenuation factor in
$z$
as well as a different force prefactor than what is found in Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989, Reference Shiokawa, Matsui and Ueda1990). In those works, the authors find that
$K_z=\alpha _1 k_{{s},{i}}$
, as opposed to our model which yields
$K_z\approx {k_{{s},{r}}\varepsilon _{{s}}/\alpha _1}=k_{{s},{i}}/\alpha _1$
. From the physical point of view, one expects larger viscous dissipation at a given
$z$
(distance from the substrate, into the liquid), when the Rayleigh angle,
$\theta _R$
, increases and tends to
$\pi /2$
(see figure 5, where the angle is measured clockwise from the normal to the substrate). This is because the leaky wave travels a longer distance along the ray (for given
$z$
) for larger
$\theta _R$
, and thus its amplitude is more attenuated. For example, if
$c \rightarrow c_{{s}}$
, we have
$\theta _R \rightarrow \pi /2$
and
$\alpha _1=1/\tan \theta _R \rightarrow 0$
(from Snell’s law, we have
$\sin (\theta _R)=c/c_{{s}} \lt 1$
). Since
$K_z$
must increase for
$\theta _R \rightarrow \pi /2$
, we have
$\tan \theta _R \rightarrow \infty$
, and then one expects
$K_z \propto 1/\alpha _1 \rightarrow \infty$
, and not
$K_z \propto \alpha _1 \rightarrow 0$
as in Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1989, Reference Shiokawa, Matsui and Ueda1990).
Appendix C. Computational implementation
Equation (3.29) is written in a form convenient for the use of COMSOL™Multiphysics PDE Coefficients Form. This package solves, by finite elements, a vectorial equation for the unknown vector
$\boldsymbol{u}=(u_1,u_2,\ldots ,u_N)^\top$
. The equation is of the form
\begin{equation} \boldsymbol{e} \frac {\partial ^2 \boldsymbol{u}}{\partial t^2}+\boldsymbol{d} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( - \boldsymbol{c} \boldsymbol{\nabla }\boldsymbol{u} - {\alpha } \, \boldsymbol{u} + {\gamma }\right ) + {\beta } \, \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{a} \, \boldsymbol{u} = \boldsymbol{f}, \end{equation}
where the coefficients of the
$N$
scalar equations are in the matrices
$\boldsymbol{e}$
,
$\boldsymbol{d}$
,
$\gamma$
,
$\boldsymbol{a}$
(of dimensions
$N\times N$
),
$\alpha$
,
$\beta$
(of dimensions
$N\times N\times n$
),
$\boldsymbol{c}$
(of dimensions
$N\times N\times n \times n$
) and the vector
$\boldsymbol{f}$
(of dimension
$N$
), where
$n$
is the spatial dimension of the problem (
$n=1,2,3$
). In index notation, this equation reads as
\begin{equation} e_{\textit{ij}}\frac {\partial ^2 u_{\!j}}{\partial \tilde {t}^2}+d_{\textit{ij}}\frac {\partial u_{\!j}}{\partial \tilde {t}}+\frac {\partial }{\partial \tilde {x}_l} \left ( -c_{\textit{ijkl}} \frac {\partial u_{\!j}}{\partial \tilde {x}_k} -\alpha _{\textit{ijl}} u_{\!j} + \gamma _{il} \right )+ \beta _{\textit{ijl}} \frac {\partial u_{\!j}}{\partial \tilde {x}_l} + a_{\textit{ij}} u_{\!j} = f_i, \end{equation}
where
$i,j=1.\ldots ,N$
and
$k,l=1,\ldots ,n$
.
The considered system, given by (3.29), contains two components
$u=(\tilde {h}, \tilde {\mathcal P})$
(
$N=2$
) and we use two equations, namely, (3.29) and (3.30) for (
$n=1)$
, corresponding to the solution depending on a single spatial variable,
$\tilde {x}$
.
In the following, we list the non-vanishing coefficients (we omit the indices
$k$
and
$l$
for brevity and consider
$\tilde {x}_1\equiv \tilde {x}$
since
$k=l=1$
):
-
(i) Row
$1$
(
$i=1$
) for (3.29)(C3)
\begin{equation} d_{11}=1,\quad c_{12}= \tilde {h}^3, \quad \gamma _1 =- \tilde {\mathcal C} \frac {3 {\mathcal S}}{8 \tilde {K}_z^4} \psi (\tilde {x}, \tilde {h}) \Big (2 \tilde {K}_z^2 \tilde {h}^2-1+ e^{2 \tilde {K}_z \tilde {h}}(1-2 \tilde {K}_z \tilde {h}) \Big ) \bigg ] \end{equation}.
-
(ii) Row
$2$
(
$i=2$
) for (3.30)(C4)
\begin{equation} c_{21} = -1, \quad a_{21}=-{\textit{Bo}},\quad a_{22}=1, \quad f_2=\frac {{\mathcal S} \tilde {C}_z}{2 \tilde {K}_z} \psi (\tilde {x}, \tilde {h}) \end{equation}.
At the domain ends, we apply Dirichlet boundary conditions
$\tilde {h}= \tilde {h}_{{p}}$
, where
$\tilde {h}_{{p}}$
is the precursor film thickness, and
$\partial \tilde {h}/\partial \tilde {x}=0$
. In order to achieve convergent numerical simulations, we discretise the domain using
$\Delta \tilde {x} =\tilde {h}_{{p}}$
.
Open access
