1. Introduction
The phenomenon of patterned wave formation on the surface of a liquid film vibrated on a surface is a classical problem, which was first studied by Faraday (Reference Faraday1831). The unstable behaviour of the film, which causes the temporal growth of surface deformations, later came to be known as the Faraday instability. The patterned waves were recognised by Lord Rayleigh as having a frequency that could be related to the forcing frequency of the substrate (Rayleigh Reference Rayleigh1883). This later came to be described by the Mathieu equation. This phenomenon finds applications in areas such as surface wave pattern formation and vibration-induced droplet ejection, especially as it relates to ultrasonic atomisation.
The present paper studies the Faraday instability of a viscous liquid film for its relevance for spray formation by a type of ultrasonic atomiser with a vibrating horn. This is an important technique for many medical and industrial applications, such as air humidification and etching of semiconductor materials. The liquid atomisation process is based on the unstable behaviour of a liquid film on a solid surface vibrating at ultrasonic frequencies. The process has attracted interest for decades, starting from the theoretical study of a vertically vibrating vessel containing an inviscid liquid with a free surface by Benjamin & Ursell (Reference Benjamin and Ursell1954). For an ideal fluid, the behaviour of the liquid–gas interface deformation in time is governed by the Mathieu equation. The solutions are presented in the form of a stability chart, depicting the surface behaviour as a function of wavenumber and vibration frequency. Benjamin & Ursell (Reference Benjamin and Ursell1954) were able to relate the forcing frequency with the wave pattern frequency for an ideal fluid with a free surface. As implied by the ideal fluid assumption, this analysis neglects the damping effect due to the fluid viscosity. The response frequency of the surface waves can be harmonic or subharmonic. A harmonic response is defined as a surface wave having the same frequency of oscillation as the vibrating plate beneath, or any integer multiple of the input frequency. A subharmonic response occurs when the surface waves oscillate at a fraction of the input frequency, typically half (or an integer fraction) of the plate’s vibration frequency.
Independently and shortly thereafter, Eisenmenger (Reference Eisenmenger1959) introduced a damping term into the Mathieu equation, known as the damped Mathieu equation. Peskin & Raco (Reference Peskin and Raco1963) used the original Mathieu equation to carry out a linear stability analysis for the ultrasonic atomisation of an inviscid liquid film to reveal the relationships between the drop size, forcing frequency, forcing amplitude, liquid-film height and physicochemical liquid properties. The influence from viscosity was included, using the results by Eisenmenger (Reference Eisenmenger1959). The expected drop size was quantified as a function of the film thickness. As it pertains to ultrasonic drop formation, the drop size is governed by the ratio of the liquid-film thickness to the vibration amplitude
$a_0$
of the substrate and by the group
$\omega ^2 a_0^3 \rho / 2 \sigma$
corresponding to a Weber number, with the angular frequency of the vibration
$\omega$
, the liquid density
$\rho$
and the surface tension
$\sigma$
. In his visualisation study, Topp (Reference Topp1973) investigated the mechanism of disintegration of a vibrating film. The study provided evidence of the wave pattern on the liquid–gas interface caused by the vibrations. Jacqmin & Duval (Reference Jacqmin and Duval1988) studied instabilities due to oscillating accelerations normal to a viscous fluid–fluid interface. Their work revealed the wavenumbers of fastest-growing waves. The wavelengths of the most unstable waves are shown to be insensitive to the differences in density and viscosity between the two fluids. Sindayihebura & Bolle (Reference Sindayihebura and Bolle1993) studied the stability behaviour of inviscid films under vibration, and Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) the corresponding viscous case. Stability charts shown depict regimes of stable and unstable film behaviour as functions of vibration frequency and wavenumber.
The work of Miles (Reference Miles1993) provided a theoretical study of Faraday waves in the context of pattern selection and provided an advanced solution for unstable capillary–gravity waves under the assumption that lateral boundary effects are negligible. That work is complementary to the comprehensive theoretical study of the Faraday instability for non-ideal fluids developed by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) to find the threshold of the stability for a single-forcing-frequency, two-fluid system. Floquet theory was applied in their solution to solve the stability problem of the linearised Navier–Stokes equations, thus accounting for viscous effects. A key restriction to their solution, which allowed them to solve for the real, positive eigenvalues of the forcing amplitude, was the selection of the response frequency of the surface waves. By setting the real part of the complex rate factor to 0, corresponding to the instability threshold, and the imaginary part of the complex rate factor to be 0 or
$\omega / 2$
, they identified the stability threshold for these cases constituting harmonic and subharmonic responses, respectively. In the stability boundary plots of Benjamin & Ursell (Reference Benjamin and Ursell1954), Kumar & Tuckerman (Reference Kumar and Tuckerman1994) and Kumar (Reference Kumar1996), the growth rates of the film surface deformation amplitude inside the regimes of unstable film behaviour are not given. The linear stability analysis performed by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) provided a basis for the theoretical analysis of the Faraday instability, to which many subsequent works have referred. Batson, Zoueshtiagh & Narayanan (Reference Batson, Zoueshtiagh and Narayanan2013) studied both experimentally and theoretically the Faraday instability threshold in cylinders at low frequencies, addressing the sidewall non-ideality by proper selection of the test liquids regarding their wetting of the container sidewall. This feature allowed those authors to detect the close connection of the experimental findings with viscous linear stability theory. Deviations seen are explained by harmonic meniscus waves. Zhao, Dietzel & Hardt (Reference Zhao, Dietzel and Hardt2019) studied the Faraday instability of a liquid layer on a lubrication film and found that the time-periodic excitation of the system creates a steady-state deformation of the bottom layer. The linear analysis was found to predict excellently the experimental findings.
In view of the stability behaviour of a vibrating liquid film, the geometrical properties of the deformation patterns formed on the liquid film surface are important. Predicting the selected planforms requires nonlinear theory. Chen & Viñals (Reference Chen and Viñals1999) formulated the standing-wave amplitude equation in a gradient form involving a Lyapunov function. The preferred film surface pattern is determined by minimisation of the latter function. Westra et al. (Reference Westra, Binks and van de Water2003) carried out experiments with a very large aspect ratio apparatus, showing the onset of various patterns with different symmetries. Good agreement was found with the predictions of the nonlinear theory by Chen & Viñals (Reference Chen and Viñals1999). Skeldon & Guidoboni (Reference Skeldon and Guidoboni2007) carried out a weakly nonlinear analysis to describe pattern selection for Faraday waves in incompressible viscous fluids for excitation with two frequencies. Results are presented in the limit of infinite liquid depth. Skeldon & Rucklidge (Reference Skeldon and Rucklidge2015) showed that weakly nonlinear theory is a useful tool to understand the transitions between different surface deformation patterns near their onset. The work by Picardo & Narayanan (Reference Picardo and Narayanan2017) showed that systems exhibiting dispersion curves with two peaks of comparable heights may be dominated by either one of the peaks, even if the linear growth rates are identical. In the present paper we base our analysis on a prescribed square deformation pattern, as seen in experiments.
Beyond studies of the Faraday instability and the phenomenon of patterned surface waves, additional literature exists on the process of film surface deformation at amplitudes large enough for atomisation to occur (Goodridge et al. Reference Goodridge, Shi, Hentschel and Lathrop1997; Goodridge Reference Goodridge1998; Yule & Al-Suleimani Reference Yule and Al-Suleimani2000; James et al. Reference James, Vukasinovic, Smith and Glezer2003; Li & Umemura Reference Li and Umemura2014). This work on vibration-induced droplet ejection has been experimental (Goodridge et al. Reference Goodridge, Shi, Hentschel and Lathrop1997; Goodridge Reference Goodridge1998; Yule & Al-Suleimani Reference Yule and Al-Suleimani2000) and numerical (James et al. Reference James, Vukasinovic, Smith and Glezer2003; Li & Umemura Reference Li and Umemura2014), rather than theoretically predicting the interfacial dynamics. A comprehensive review of capillary effects on surface waves is due to Perlin & Schultz (Reference Perlin and Schultz2000).
In the present paper we study the linear stability behaviour of a viscous liquid film on a plane horizontal vibrating substrate. The relevance of the study is for ultrasonic atomisation, as reviewed above. Though the results presented are focused on the range of frequencies associated with ultrasonic atomisation, the present method is general to a liquid film vibrated on a flat plate across scales. The aim is to quantify the growth rate and the oscillation frequency of the surface deformations of the liquid film as functions of all the influencing parameters. In the following section, the problem is formulated, and the governing equations with their boundary conditions are derived. Section 3 presents the solutions for the disturbance velocity components and pressure. The dispersion relation of the film is derived from these solutions in § 4. In § 5, the solutions are evaluated, and the results are discussed. The paper ends with the conclusions in § 6.
2. Formulation of the problem: governing equations
We study the linear stability of a liquid film with a free surface in contact with an ambient inert gaseous medium. The film is placed on a vibrating solid substrate used for destabilising the film and for the formation of a fine spray. The corresponding technical process is known as ultrasonic atomisation of a liquid. The aim of the present work is to provide a linear stability analysis of the liquid film. In the literature, there is a paper by Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) that presents such an analysis. Verification of the calculations, however, shows that an important part of the solution of a differential equation was missed.
In the present analysis, the liquid film with thickness
$h$
is treated as incompressible and viscous Newtonian. The equations of change are formulated in a Cartesian coordinate system moving with the vibrating substrate surface. To do this, we account for the sinusoidal vibrating motion of the substrate with amplitude
$a_0$
and angular frequency
$\omega$
by its acceleration
$a_0 \omega ^2 \cos \omega t$
adding to the acceleration due to body forces. The term formed by the transformation from the laboratory-fixed to the co-moving coordinate system is then merged with the pressure gradient. The equations of motion are non-dimensionalised by
$h$
,
$(h/g)^{1/2}$
,
$(gh)^{1/2}$
and
$\rho g h$
for length, time, velocity and pressure, respectively. The vector of gravitational acceleration points in the negative
$z$
direction, normal to the undisturbed free liquid-film surface. In the moving reference frame, the undisturbed liquid film is at rest, i.e. all the liquid velocities we obtain as a result of our analysis are due to the disturbance by the shaking. The system is sketched in figure 1. With the disturbance velocity vector
$(u,v,w)$
, the non-dimensional continuity equation reads
\begin{equation} \frac {\partial u}{\partial x} + \frac {\partial v}{\partial y} + \frac {\partial w}{\partial z} = 0. \end{equation}
Assuming small liquid-film velocity disturbances, the momentum equations are linearised. The non-dimensional linearised balance equations for the
$x$
,
$y$
and
$z$
momentum read
\begin{align} \frac {\partial u}{\partial t} & = - \frac {\partial p}{\partial x} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} + \frac {\partial ^2 u}{\partial z^2} \right ]\!, \\[-9pt] \nonumber \end{align}
\begin{align} \frac {\partial v}{\partial t} & = - \frac {\partial p}{\partial y} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 v}{\partial x^2} + \frac {\partial ^2 v}{\partial y^2} + \frac {\partial ^2 v}{\partial z^2} \right ]\!, \\[-9pt] \nonumber \end{align}
\begin{align} \frac {\partial w}{\partial t} & = - \frac {\partial p}{\partial z} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 w}{\partial x^2} + \frac {\partial ^2 w}{\partial y^2} + \frac {\partial ^2 w}{\partial z^2} \right ] - \left ( 1 - \textit{Fr} \cos \varOmega t \right )\!, \\[9pt] \nonumber \end{align}
respectively, where
\begin{equation} \textit{Re} = \frac {(gh)^{1/2} h}{\nu }, \quad \textit{Fr} = \frac {a_0 \omega ^2}{g}, \quad \varOmega = \omega \left ( \frac {h}{g} \right )^{1/2} \end{equation}
are the Reynolds number, the Froude number and the non-dimensionalised angular frequency of the substrate motion.
Sketch of the geometry of the film with surface deformations on the substrate (adapted from Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998)). The term
$a_0 \omega ^2 \cos \omega t$
denotes the acceleration of the system due to the vibrations.
Defining a modified pressure as
\begin{equation} p' = p + \left ( 1 - \textit{Fr} \cos \varOmega t \right ) z, \end{equation}
we rewrite the set of momentum equations as
\begin{align} \frac {\partial u}{\partial t} & = - \frac {\partial p'}{\partial x} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 u}{\partial x^2} + \frac {\partial ^2 u}{\partial y^2} + \frac {\partial ^2 u}{\partial z^2} \right ]\!, \\[-9pt] \nonumber\end{align}
\begin{align} \frac {\partial v}{\partial t} & = - \frac {\partial p'}{\partial y} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 v}{\partial x^2} + \frac {\partial ^2 v}{\partial y^2} + \frac {\partial ^2 v}{\partial z^2} \right ]\!, \\[-9pt] \nonumber\end{align}
\begin{align} \frac {\partial w}{\partial t} & = - \frac {\partial p'}{\partial z} + \frac {1}{\textit{Re}} \left [ \frac {\partial ^2 w}{\partial x^2} + \frac {\partial ^2 w}{\partial y^2} + \frac {\partial ^2 w}{\partial z^2} \right ]\!. \\[9pt] \nonumber\end{align}
The solution of these equations of motion is developed for a laterally infinite liquid film, subject to the boundary conditions
\begin{equation} u = v = w = 0 \quad \mbox{at} \quad z = -1. \end{equation}
Furthermore, the condition
\begin{equation} \frac {\partial w}{\partial z} = 0 \quad \mbox{at} \quad z = -1 \end{equation}
applies, which follows from the continuity (2.1). The
$z$
position of the deformed liquid-film surface, which at equilibrium is located at
$z = 0$
, is denoted as
$\zeta (x, y, t)$
. The kinematic boundary condition at the free surface of the film states that
\begin{equation} w = \frac {\partial \zeta }{\partial t} \quad \mbox{at} \quad z = 0. \end{equation}
The ambient medium is treated as dynamically inert, so that the viscous shear stress at the film surface is zero. The linear zero-shear-stress boundary conditions read
\begin{align} \frac {\partial v}{\partial z} + \frac {\partial w}{\partial y} & = 0 \quad \mbox{at} \quad z = 0, \\[-9pt] \nonumber\end{align}
\begin{align} \frac {\partial w}{\partial x} + \frac {\partial u}{\partial z} & = 0 \quad \mbox{at} \quad z = 0. \\[9pt] \nonumber\end{align}
The normal stress on the liquid side of the liquid–gas interface must likewise be zero, since the ambient medium is dynamically inert. The influence from the surface curvature is accounted for by a capillary stress due to the surface tension
$\sigma$
of the liquid against the vacuum. The normal-stress boundary condition therefore reads
\begin{equation} -p + \frac {2}{\textit{Re}} \frac {\partial w}{\partial z} - \frac {1}{\textit{We}} \left [ \frac {\partial ^2 \zeta }{\partial x^2} + \frac {\partial ^2 \zeta }{\partial y^2} \right ] = 0 \quad \mbox{at} \quad z = 0, \end{equation}
where the Weber number is defined as
$\textit{We} = g h^2 \rho / \sigma$
. Yet, the condition
\begin{equation} \frac {\partial p'}{\partial z} = \frac {1}{\textit{Re}} \frac {\partial ^2 w}{\partial z^2} \quad \mbox{at} \quad z = -1 \end{equation}
states that, at the bottom of the liquid film, the modified pressure gradient in the
$z$
direction is in equilibrium with the gradient of the normal stress.
Top view of the surface of a liquid film on the vibrating front surface of an ultrasonic atomiser. The wave crests exhibit a square pattern. Vibration frequency, 20 kHz; transducer diameter, 14.5 mm. Reprinted from Ramisetty et al. (Reference Ramisetty, Pandit and Gogate2013) with permission from Elsevier.
3. Solutions of the governing equations
In the present section, solutions of the governing equations are determined, so as to describe the three disturbance velocity components in space,
$u,\,v$
and
$w$
, and the disturbance pressure
$p'$
. While the present analysis describes a three-dimensional system such as is seen in figure 2, restricting the solution to two-dimensional Cartesian coordinates would also be possible by reformulating the
$x$
and
$y$
dependencies in the pressure and velocity fields. Ultimately this would not change the amplitude functions or the resulting dispersion relation.
3.1. The disturbance pressure
Forming the divergence of the vectorial momentum equation, and using the continuity equation, shows that the pressure
$p'$
satisfies a Laplace equation. In view of the normal-stress boundary condition, the pressure must match the shape of the deformed liquid surface. Shapes of deformed liquid-film surfaces found in the experiment may look like that in figure 2, showing a square-wave pattern, despite the axially symmetric form of the vibrating solid transducer (Ramisetty, Pandit & Gogate Reference Ramisetty, Pandit and Gogate2013). Similar square-wave patterns were shown by Sindayihebura et al. (Reference Sindayihebura, Dobre and Bolle1997) for the vibration frequencies
$32$
and
$50$
kHz in figure 1 of their paper also. It is important to note that the present analysis relies on linear stability theory, assuming a laterally infinite domain. This corresponds to the process we are looking at, which is represented by the photograph in figure 2 (Ramisetty et al. Reference Ramisetty, Pandit and Gogate2013). The liquid film spreads on the front face of a vibrating horn. That horn surface has no sidewalls, which is why our analysis does not account for them. This approach therefore neglects the influence of sidewalls and spatial confinement, which are present in other experimental set-ups and play a crucial role in determining the final planform pattern. Linear theory alone cannot predict the experimentally observed square planforms, which arise from nonlinear interactions by higher-order effects. Nonlinear studies, such as that by Skeldon & Rucklidge (Reference Skeldon and Rucklidge2015), offer important insight into the mechanisms governing pattern formation in confined systems. This is outside the scope of the present linear analysis. Such wave patterns are consistent with a mathematical description of the surface deformation
$\zeta$
against the equilibrium position at
$z=0$
by
$\zeta \propto ( \cos kx + \cos ky )$
, where
$k = 2 \pi h / \lambda$
is the non-dimensional wavenumber, which, in view of the experimental findings, has the same value in both coordinate directions
$x$
and
$y$
(Rayleigh Reference Rayleigh1883). A solution of the Laplace equation for the pressure disturbance
$p'(x,y,z,t)$
is therefore formulated as
\begin{equation} p' (x, y, z, t) = \left ( \cos kx + \cos ky \right ) \big ( C_1 (t) {\rm e}^{kz} + C_2 (t) {\rm e}^{-kz} \big ), \end{equation}
from where we derive the definition of the amplitude function
${\hat p} (z,t)$
of
$p'$
as per
\begin{equation} p' (x, y, z, t) = \left ( \cos kx + \cos ky \right ) {\hat p} (z,t). \end{equation}
Consequently, a solution
$w (x,y,z,t)$
of the
$z$
momentum equation (2.9) and a corresponding shape
$\zeta (x,y,t)$
of the film surface must be of the forms
\begin{align} w (x, y, z, t) & = \left ( \cos kx + \cos ky \right ) \,{\hat h} (z,t), \\[-9pt] \nonumber\end{align}
\begin{align} \zeta (x, y, t) & = \left ( \cos kx + \cos ky \right ) \,{\hat \zeta } (t) \\[9pt] \nonumber\end{align}
in order to satisfy the momentum equations and the boundary conditions. The amplitude functions
${\hat p}(z,t)$
,
${\hat h} (z,t)$
and
${\hat \zeta } (t)$
, governing the
$z$
and time dependencies of the solutions, remain to be determined.
The time-dependent coefficients
$C_1 (t)$
and
$C_2 (t)$
in the solution (3.1) are determined by the conditions (2.15) and (2.16) and read
\begin{align} C_1 (t) & = \big ( k^2 / \textit{We} + 1 - \textit{Fr} \cos \varOmega t \big ) {\hat \zeta } (t) \frac {{\rm e}^k}{2 \cosh k} + \frac {{\rm e}^k}{\textit{Re} \cosh k} \left . \left ( \frac {\partial {\hat h}}{\partial z} \right ) \right |_{z=0}\nonumber\\& \quad + \frac {1}{2 \,\textit{Re}\, k \cosh k} \left . \left ( \frac {\partial ^2 {\hat h}}{\partial z^2} \right ) \right |_{z=-1}, \\[-9pt] \nonumber\end{align}
\begin{align} C_2 (t) & = \big ( k^2 / \textit{We} + 1 - \textit{Fr} \cos \varOmega t \big ) {\hat \zeta } (t) \frac {{\rm e}^{-k}}{2 \cosh k} + \frac {{\rm e}^{-k}}{\textit{Re} \cosh k} \left . \left ( \frac {\partial {\hat h}}{\partial z} \right ) \right |_{z=0} \nonumber \\& \quad - \frac {1}{2\, \textit{Re}\, k \cosh k} \left . \left ( \frac {\partial ^2 {\hat h}}{\partial z^2} \right ) \right |_{z=-1}. \\[9pt] \nonumber\end{align}
The solution for the pressure
$p'$
follows as
\begin{eqnarray} p' (x, y, z, \tau ) & = & \left \{ \left [ \big ( k^2/\textit{We} + 1 - \textit{Fr} \cos 2 \tau \big ) {\hat \zeta } (\tau ) + \frac {2}{\textit{Re}} \left . \left ( \frac {\partial {\hat h}}{\partial z} \right ) \right |_{z=0} \right ] \frac {\cosh k(z+1)}{\cosh k} \right . \nonumber \\& + & \quad \left . \frac {1}{k \textit{Re}} \left . \left ( \frac {\partial ^2 {\hat h}}{\partial z^2} \right ) \right |_{z=-1} \frac {\sinh kz}{\cosh k} \right \} ( \cos kx + \cos ky ) , \end{eqnarray}
where we have defined
$\tau = \varOmega t / 2$
. This solution (3.7) determines the amplitude function
${\hat p}(z,t)$
as the term in the curled brackets. The function
$\hat {h}(z,t)$
, which appears in the equation for
$\hat {p}(z,t)$
, is yet to be determined. In the following section, expressions for the amplitude functions of the velocity components are derived.
3.2. The disturbance velocity components
As a next step of the analysis, the
$z$
momentum equation (2.9) is used to determine the disturbance velocity
$w$
in the
$z$
direction. For this purpose, the functions
${\hat h} (z,\tau )$
and
${\hat \zeta } (\tau )$
are formulated such that their time dependencies are expressed as Fourier series of
$\pi$
-periodic functions multiplied by an exponential function exhibiting a complex rate factor consisting of an angular frequency and a rate of growth or decay in time. The
$\pi$
periodicity corresponds to the function
$\cos 2 \tau$
in the solution (3.7). We formulate
\begin{align} {\hat h} (z, \tau ) & = \sum \limits _{i = - \infty }^{\infty } H_i (z) {\rm e}^{j m_i \tau } \quad \mbox{and} \\[-9pt] \nonumber \end{align}
\begin{align} {\hat \zeta } (\tau ) & = \sum \limits _{i = - \infty }^{\infty } N_i {\rm e}^{j m_i \tau }, \\[9pt] \nonumber \end{align}
where we denote the imaginary unit as
$j = \sqrt {-1}$
and formulate the rate factor
\begin{equation} { m_i = 2 i - j \alpha } \end{equation}
to account for the
$\pi$
periodicity of the solutions and to introduce the rate factor
$\alpha$
representing damping or growth in time of the disturbances caused by the substrate vibrations. Positive values of the real part of
$\alpha$
indicate growth of the disturbance in time. Together with this representation of the time dependency by exponential functions, we write
$\cos 2 \tau$
as
\begin{equation} \cos 2 \tau = \frac {1}{2} \big ( {\rm e}^{j 2 \tau } + {\rm e}^{-j 2 \tau } \big ). \end{equation}
With the formulations (3.2) and (3.3), the
$z$
momentum equation (2.9) turns into the form
\begin{equation} \frac {\partial {\hat h}}{\partial \tau } \frac {\varOmega }{2} = - \frac {\partial {\hat p}}{\partial z} + \frac {1}{\textit{Re}} \left [ -k^2 {\hat h} + \frac {\partial ^2 {\hat h}}{\partial z^2} \right ]\!, \end{equation}
where the functions
$(\cos kx + \cos ky)$
were cancelled. Using the amplitude function
$\hat p$
from the solution (3.7) and the formulations (3.8) and (3.9), the equation becomes
\begin{eqnarray} & & \sum \limits _{i = - \infty }^{\infty } \left [ H^{\prime\prime}_i(z) - \left ( k^2+jm_i \frac {\varOmega \textit{Re}}{2} \right ) H_i (z) \right ] {\rm e}^{j m_i \tau } \nonumber \\& & =\left [ 2k \sum \limits _{i = - \infty }^{\infty } H^{\prime}_i (0) {\rm e}^{j m_i \tau } + \left ( k\,\textit{Re} \left ( \frac {k^2}{\textit{We}} + 1 \right ) \right . \right . \nonumber \\& & \left . \left . \quad - \frac {k\,\textit{Re}}{2} \textit{Fr} \big({\rm e}^{j 2 \tau } + {\rm e}^{-j 2 \tau } \big) \right ) \sum \limits _{i = - \infty }^{\infty } N_i {\rm e}^{j m_i \tau } \right ] \frac {\sinh k(z+1)}{\cosh k} \nonumber \\& & \quad + \sum \limits _{i = - \infty }^{\infty } H^{\prime\prime}_i(-1) {\rm e}^{j m_i \tau } \frac {\cosh kz}{\cosh k}. \end{eqnarray}
We rewrite the term with
$\cos 2 \tau$
in the third line of this equation as
\begin{eqnarray} \big({\rm e}^{j 2 \tau } + {\rm e}^{-j 2 \tau } \big) \sum \limits _{i = - \infty }^{\infty } N_i {\rm e}^{j m_i \tau } & = & \sum \limits _{i = - \infty }^{\infty } N_i {\rm e}^{(j2(i+1)+\alpha )\tau } + \sum \limits _{i = - \infty }^{\infty } N_i {\rm e}^{(j2(i-1)+\alpha )\tau } \nonumber \\& = & \sum \limits _{i = - \infty }^{\infty } N_{i-1} {\rm e}^{(j2i+\alpha )\tau } + \sum \limits _{i = - \infty }^{\infty } N_{i+1} {\rm e}^{(j2i+\alpha )\tau } \nonumber \\& = & \sum \limits _{i = - \infty }^{\infty } \big ( N_{i-1} + N_{i+1} \big ) {\rm e}^{j m_i \tau }, \end{eqnarray}
where the second equality is due to summation index transforms in the two sums. This allows the differential equation (3.13) to be represented by the following equation for each term in the sums:
\begin{eqnarray} & & H^{\prime\prime}_i(z) - \xi _i^2 H_i (z) = \left [ 2k H^{\prime}_i (0) + k\,\textit{Re} \left ( \frac {k^2}{\textit{We}} + 1 \right ) N_i \right . \nonumber \\& - & \left . \frac {k\,\textit{Re}}{2} \textit{Fr} \big ( N_{i-1} + N_{i+1} \big ) \right ] \frac {\sinh k(z+1)}{\cosh k} + H^{\prime\prime}_i(-1) \frac {\cosh kz}{\cosh k}, \end{eqnarray}
where we have defined
\begin{equation} {\xi _i^2 = k^2 + j m_i \varOmega \textit{Re}/2 \equiv k^2 + j m_i \omega h^2 / (2 \nu )}. \end{equation}
This differential equation has the form
\begin{equation} H^{\prime\prime}_i(z) - \xi _i^2 H_i (z) = \beta \frac {\sinh k(z+1)}{\cosh k} + \gamma \frac {\cosh kz}{\cosh k}, \end{equation}
where the coefficients
$\beta$
and
$\gamma$
represent the factors in front of the fractions of hyperbolic functions in (3.15). The corresponding homogeneous differential equation has the general solution
\begin{equation} H_{i,h}(z) = D_1 {\rm e}^{\xi _i z} + D_2 {\rm e}^{- \xi _i z}, \end{equation}
with the complex coefficients
$D_1$
and
$D_2$
. A particular solution of the inhomogeneous differential equation is found as
\begin{equation} H_{i,p}(z) = \frac {\beta }{k^2 - \xi _i^2} \frac {\sinh k(z+1)}{\cosh k} + \frac {\gamma }{k^2 - \xi _i^2} \frac {\cosh kz}{\cosh k}. \end{equation}
The general solution of the differential equation (3.15) reads
\begin{equation} H_i (z) = H_{i,h}(z) + H_{i,p}(z). \end{equation}
The boundary conditions (2.10) and (2.11) state that both the
$z$
component of the velocity and the
$z$
component of its gradient must be zero at the bottom of the liquid film,
$z = -1$
. These conditions allow the two coefficients
$D_1$
and
$D_2$
in the solution of the homogeneous equation to be determined. We obtain
\begin{equation} D_1 = \left ( A_1 + A_2 \right ) \frac {{\rm e}^{\xi _i}}{2}, \quad D_2 = \left ( A_1 - A_2 \right ) \frac {{\rm e}^{-\xi _i}}{2}, \end{equation}
where
$A_1$
and
$A_2$
are given as
\begin{equation} A_1 = - \frac {\gamma }{k^2 - \xi _i^2}, \quad A_2 = \frac {k}{\xi _i \big( k^2 - \xi _i^2 \big) \cosh k} ( \gamma \sinh k - \beta ). \end{equation}
We therefore rewrite the solution for
$H_i (z)$
as
\begin{align} H_i (z) & = A_1 \cosh \xi _i(z+1) + A_2 \sinh \xi _i(z+1) + \frac {\beta }{k^2 - \xi _i^2} \frac {\sinh k(z+1)}{\cosh k} \nonumber \\ & \quad + \frac {\gamma }{k^2 - \xi _i^2} \frac {\cosh kz}{\cosh k}. \end{align}
Using this solution for
$H_i$
allows an equation for its first derivative at
$z = 0$
to be determined. The second derivative at
$z = -1$
, however, remains undetermined. Comparing the solution with (39) from Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) reveals that the latter paper represented the surface wave pattern by a product of a Bessel function
$J_n(kr)$
and a cosine function
$\cos n \theta$
, with constant wavenumbers
$k$
and
$n$
. The square surface wave pattern seen in the experiment is compared in figure 3 with that predicted by Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) for the value of
$n=5$
as an example. The figure shows that the solution based on cylindrical functions and a cosine does not represent the pattern that is experimentally observed in cases of resonance where unstable behaviour of the film is expected. The square pattern is known as characteristic for a resonant motion of the vibrating substrate. The pattern exhibiting a concentric circular shape is formed in the experiment for states off resonance, but it does not exhibit a wavy form in the azimuthal direction (Sindayihebura Reference Sindayihebura1995).
Three-dimensional view of a vibrating film surface displaying (a) a generic square standing-wave pattern, which qualitatively mimics what is seen in figure 2, and (b) a pattern of the form
$\cos n \theta \, J_n (kr)$
, with the polar angle
$\theta = \arctan \,y/x$
, as predicted by Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998), the latter for the radial wavenumber
$k = 1.1318$
, derived from the first zero of the Bessel function
$J_5(kr)$
, and the angular wavenumber
$n=5$
as an example. The coordinate axes in the horizontal plane display
$x$
and
$y$
.
For further developing the analysis, we formulate the disturbance velocity component in the
$x$
direction,
$u(x, y, z, t)$
. For formulating its dependencies on the
$x$
and
$y$
coordinates, we use the zero-shear-stress boundary condition (2.14). From this condition it follows that
$u$
must depend on
$x$
as per
$\sin kx$
, while it does not depend on the coordinate
$y$
. We therefore formulate
\begin{equation} u (x, y, z, t) = \sin kx \,{\hat g} (z,\tau ). \end{equation}
Introducing this formulation into the
$x$
momentum equation (2.7) yields the differential equation
\begin{equation} \frac {\partial {\hat g}}{\partial \tau } \frac {\varOmega }{2} = k {\hat p} + \frac {1}{\textit{Re}} \left ( - k^2 {\hat g} + \frac {\partial ^2 {\hat g}}{\partial z^2} \right ) \end{equation}
for the amplitude function
${\hat g} (z,\tau )$
, where the function
$\sin kx$
was cancelled. The separate
$z$
-coordinate and time dependencies of the function
${\hat g} (z,\tau )$
are formulated, analogous to (3.8), as
\begin{equation} {\hat g} (z, \tau ) = \sum \limits _{i = - \infty }^{\infty } G_i (z) {\rm e}^{j m_i \tau }. \end{equation}
The so formulated function
${\hat g} (z,\tau )$
must satisfy (3.25). For determining
$G_i(z)$
, we follow a way analogous to that applied for getting
$H_i(z)$
. The solution for
$G_i (z)$
, formulated with the same integration constants
$A_1$
and
$A_2$
as in the solution (3.23) for
$H_i (z)$
, reads
The so formulated disturbance velocity component
$u$
satisfies the boundary conditions (2.10) and (2.14).
By virtue of the boundary condition (2.13), the disturbance velocity component
$v(x, y, z, t)$
in the direction of the coordinate
$y$
is found to depend on the
$y$
coordinate as per
$\sin ky$
. The velocity component is then searched as a solution of the
$y$
momentum balance (2.8). The equation shows that the solution for the function governing the
$z$
and
$\tau$
dependencies of
$v$
is the same as for the disturbance velocity component
$u$
. The formulation of the solution for
$v$
therefore reads
\begin{equation} v (x, y, z, t) = \sin ky \,{\hat g} (z,\tau ). \end{equation}
The so formulated disturbance velocity component
$v(x, y, z, t)$
satisfies the boundary conditions (2.10) and (2.13).
The above solutions for the disturbance velocity components
$u,\,v$
and
$w$
involve the as yet undetermined quantity
$\gamma \equiv H^{\prime\prime}_i(-1)$
. This quantity is determined, requiring that the velocity components must satisfy the continuity equation (2.1). Substituting the solutions into this equation yields the requirement that
\begin{equation} - \gamma + \frac {k}{\xi _i} \frac {\tanh \xi _i}{\cosh k} \left ( \gamma \sinh k - \beta \right ) + \frac {2 k^2}{\big(k^2 + \xi _i^2 \big) \cosh k \cosh \xi _i} \left ( \beta \sinh k + \gamma \right ) = 0. \end{equation}
Furthermore, another equation involving
$\beta$
is found by formulating
$\beta = 2 k H^{\prime}_i (0) + J$
, substituting the expression for
$H^{\prime}_i(0)$
using (3.23) and solving for
$\beta$
, where we have abbreviated
\begin{equation} J := k \textit{Re} \left ( \frac {k^2}{\textit{We}} + 1 \right ) N_i - \frac {k \textit{Re}}{2} \textit{Fr} \big ( N_{i-1} + N_{i+1} \big ). \end{equation}
This second relation involving
$\beta$
reads
\begin{equation} \beta = \frac {(k \tanh k \cosh \xi _i - \xi _i \sinh \xi _i) \gamma + J \big(k^2 - \xi _i^2 \big) / 2k}{k (\cosh \xi _i / \cosh k - 1) + \big(k^2 - \xi _i^2 \big) / 2k}. \end{equation}
Equating the expressions for
$\beta$
emerging from (3.29) and (3.31) allows the quantity
$\gamma \equiv H^{\prime\prime}_i(-1)$
to be determined as
\begin{align} \gamma & = \frac {k \big(k^2 - \xi _i^2 \big) \big[ \big(k^2+\xi _i^2 \big) \sinh \xi _i - 2 k \xi _i \sinh k \big]\,J}{ \xi _i \big(5 k^4 + 2 k^2 \xi _i^2 + \xi _i^4 \big)\! \cosh k \cosh \xi _i - k \big[4 k \xi _i \big(k^2 + \xi _i^2\big) + \big(k^4 + 6 k^2 \xi _i^2 + \xi _i^4\big)\! \sinh k \sinh \xi _i \big] } \nonumber \\ & =: F\,J. \end{align}
For the first-order derivative of
$H_i(z)$
at
$z = 0$
we obtain
\begin{eqnarray} H^{\prime}_i(0) & = & \frac {(k \cosh \xi _i \tanh k - \xi _i \sinh \xi _i) F - (\cosh \xi _i / \cosh k - 1) k }{\big(k^2 - \xi _i^2 \big)+ 2 k^2 (\cosh \xi _i / \cosh k - 1)} J \nonumber \\ & =: & B\,J. \end{eqnarray}
With these solutions, the three disturbance velocities and the disturbance pressure satisfying the boundary conditions are known.
4. Stability analysis
In the present section, the boundary condition (2.12) at
$z=0$
is used with (3.4) and (3.9) to determine the stability behaviour of the liquid film in response to the vibrational excitation by the periodic substrate motion. Using with this boundary condition the solutions for the
$z$
-disturbance velocity component and the formulation for the liquid-film surface deformation in time, we obtain the equation
\begin{equation} H_i(0) = j m_i \frac {\varOmega }{2} N_i. \end{equation}
Using the solution (3.23) for
$H_i(z)$
, with some algebra we arrive at the following equation determining the coefficients
$N_i$
in the formulation of the time behaviour of the film surface deformation
$\hat {\zeta } (\tau )$
:
\begin{equation} \left [ \frac {2}{\textit{Fr}} \left ( \frac {k^2}{\textit{We}} + 1 \right ) - \frac {(m_i \varOmega / 2)^2}{\frac {k}{2} \textit{Fr} F_1} \right ] N_i = N_{i-1} + N_{i+1}, \end{equation}
where we have used the abbreviation
\begin{equation} F_1 = F\,K - \frac {k}{\xi _i} \frac {\sinh \xi _i}{\cosh k} + \tanh k + 2k \left ( - \frac {k}{\xi _i} \frac {\sinh \xi _i}{\cosh k} + \tanh k \right ) B \end{equation}
with
\begin{equation} K = - \cosh \xi _i + \frac {k}{\xi _i} \sinh \xi _i \tanh k + \frac {1}{\cosh k}. \end{equation}
Denoting the factor in square brackets in front of
$N_i$
in (4.2) as
$\eta _i$
, we obtain the homogeneous system of equations
\begin{equation} N_{i-1} - \eta _i N_i + N_{i+1} = 0 \end{equation}
determining the expansion coefficients
$N_i$
in the time dependency of the liquid surface deformation amplitude
$\hat \zeta (\tau )$
. The coefficient
$\eta _i$
as per the contents of the square brackets in front of
$N_i$
in (4.2) may be simplified, since the arguments
$k$
and
$\xi _i$
of the hyperbolic functions are large,
$\approx O(50)$
. Replacing the functions
$\tanh k$
and
$\tanh \xi _i$
by unity, and neglecting one term proportional to
$1/(\cosh k \cosh \xi _i)$
yields the simple form
\begin{equation} \eta _i {\approx } \frac {2}{\textit{Fr}} \left ( \frac {k^2}{\textit{We}} + 1 \right ) + \frac {2}{\textit{Fr}} \frac {1}{k \textit{Re}^2} \big ( \xi _i^4 + 2 k^2 \xi _i^2 - 4 k^3 \xi _i + k^4 \big ) . \end{equation}
Further expanding the
$\xi _i$
in this expression, we see the explicit dependence of
$\eta _i$
on the rate factor
$\alpha$
:
\begin{align} \eta _i & \approx \frac {2}{\textit{Fr}} \left [ \frac {4k^3}{\textit{Re}^2}+ \frac {k^2}{\textit{We}}-\frac {2k^2}{\textit{Re}^2}\sqrt {4k^2+2\textit{Re} \varOmega (\alpha +j2i)}+\frac {2k \varOmega }{\textit{Re}}(\alpha +j2i)\right.\nonumber \\ & \quad \left. +\frac {\varOmega ^2}{4 k}(\alpha +j2i)^2+1 \right ]\!. \end{align}
Our coefficient
$\eta _i$
corresponds to the factor
$c_{2i}$
in (40) of Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998), but with differences due to the different derivations of the 1998 paper. Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) rewrote their equation, dividing by
$N_i$
and putting their
$c_{2i}$
on the left-hand side. Writing the equation for different values of the index
$i$
allowed them to derive expressions for their
$N_0 / N_{-2}$
and
$N_{-2} / N_0$
by successively eliminating ratios of
$N_i$
with other values of
$i$
. The results were two different expressions denoted
$G_0$
and
$H_0$
with continued fractions consisting of their
$c_{2i}$
only. The expressions
$G_0$
and
$H_0$
are inverse to each other. The requirement to these fractions was therefore that their product yields unity, i.e.
$G_0 H_0 = 1$
. This equation is the dispersion relation of the vibrating liquid film used for determining the growth or damping rate of the disturbances. In their paper, Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) do not disclose the details of their calculation of the growth rate.
In our present study we adopt the approach by Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998). Writing down the system of (4.5) for varying positive and negative values of the index
$i$
, we arrive at the two equations (Sindayihebura & Bolle Reference Sindayihebura and Bolle1998)
\begin{align} \frac {N_{-1}}{N_0} & = \eta _0 - \frac {1}{\eta _1-\frac {1}{\eta _2-\frac {1}{\eta _3 - \frac {1}{\ldots }}}}, \\[-9pt] \nonumber\end{align}
\begin{align} \frac {N_{0}}{N_{-1}} & = \eta _{-1} - \frac {1}{\eta _{-2}-\frac {1}{\eta _{-3}-\frac {1}{\eta _{-4} - \frac {1}{\ldots }}}} \\[9pt] \nonumber\end{align}
exhibiting continuous fractions on the right-hand sides. Since the quantities on the left-hand sides of these equations are inverse to each other, the product of the right-hand sides equals unity, so that the equation
\begin{equation} \left ( \eta _0 - \frac {1}{\eta _1-\frac {1}{\eta _2-\frac {1}{\eta _3 - \frac {1}{\ldots }}}} \right ) \boldsymbol{\cdot }\left ( \eta _{-1} - \frac {1}{\eta _{-2}-\frac {1}{\eta _{-3}-\frac {1}{\eta _{-4} - \frac {1}{\ldots }}}} \right ) = 1, \end{equation}
determining the rate factor
$\alpha$
, emerges as the dispersion relation of the liquid film on the vibrating substrate. The values of
$\alpha$
as a function of the wavenumber
$k$
converge to a stable shape with increasing maximum index
$i$
of
$\eta _i$
to which the continuous fractions are carried in the evaluation of the dispersion relation. The software package Wolfram Mathematica was employed to solve this dispersion relation for the complex rate factor
$\alpha$
using the ‘FindRoot’ function. The results presented in the following section were obtained by limiting the values of the summation index
$i$
to the range from a minimum value of
$-9$
to a maximum value of
$8$
.
5. Results and discussion
5.1. Properties of the dispersion relation
The dispersion relation (4.10) derived above was evaluated for quantifying the rate factor
$\alpha$
as a complex quantity
$\alpha _{\textit{real}} + j \alpha _{\textit{imag}}$
, varying all the influencing parameters. In contrast to the earlier work of Kumar & Tuckerman (Reference Kumar and Tuckerman1994), where the imaginary part of
$\alpha$
was prescribed, in our formulation it is solved for explicitly. This enables our analysis to capture Hopf bifurcations, which were not previously captured in the Kumar & Tuckerman (Reference Kumar and Tuckerman1994) framework. Hopf bifurcations describe a point where a system begins to oscillate (has a non-zero frequency) as stability is lost. As a result, complex conjugate pairs emerge as solutions in the unstable region where
$\alpha _{\textit{real}}\geqslant 0$
. This concept is highlighted in this section.
More recently, Ignatius et al. (Reference Ignatius, Dinesh, Dietze and Narayanan2024) employed Floquet analysis to examine the temporal evolution of linear perturbations under parametric forcing without prescribing values of the response frequency of the surface waves. However, like Kumar & Tuckerman (Reference Kumar and Tuckerman1994), they focused on the threshold of stability, solving for neutral conditions where the real part of the eigenvalue (the growth rate) vanishes. Our analysis, while including discussion of the neutral stability curve, also presents non-zero growth rates and the response frequency of the surface waves, which both emerge naturally as parts of the solution.
We first analyse the rate factors
$\alpha$
for liquid water and ethanol films on the vibrating surface. The liquid material data used in the analysis are listed in table 1. The real and imaginary parts of
$\alpha$
for different vibration amplitudes
$a_0$
for water films of different thicknesses are shown in figure 4. The corresponding results for ethanol films are shown in figure 5. The graphs show that the regimes of unstable film behaviour, characterised by positive
$\alpha _{\textit{real}}$
, widen with increasing
$a_0$
. This means that increasing vibration amplitude leads to the widening of the unstable region, which grows to include higher and lower wavenumbers with unstable film behaviour at larger values of
$a_0$
. The growth of the wavenumber ranges is, however, limited laterally since two closed contours do not intersect. For a given film thickness, an increase of the amplitude
$a_0$
raises the maximum film deformation growth rate. In contrast, for a given amplitude
$a_0$
, an increase of the film thickness leaves the maximum non-dimensional growth rate of the film surface deformation unaffected. The most dangerous state of vibration of the liquid film is the state (
$a_0, k$
) exhibiting the largest value of
$\alpha _{\textit{real}}$
. The results in figures 4 and 5 show that, for the values of the amplitude
$a_0$
represented there, these values always occur in the first lobe of positive
$\alpha _{\textit{real}}$
, i.e. in the unstable regime with the smallest wavenumbers. A comparison of the data in figures 4 and 5 for the different liquids furthermore shows that the disturbance growth rates for ethanol are larger than those for water, which is due to the dominant effect of surface tension, while the dampening effect of the viscosity is minor, since the dynamic viscosities of the two liquids differ by no more than
$19\,\%$
, while the surface tensions are different by a factor of
$3.3$
. For flat liquid films, high surface tension hinders the deformation, while it may enhance deformation of round liquid jets.
Material properties of water and ethanol at
$20\, ^\circ \text{C}$
.
(a,c,e) Real and (b,d,f) imaginary parts of the rate factor
$\alpha$
for (a,b)
$100$
μm and (c,d)
$500$
μm thick water films vibrating at
$20\,\rm \rm kHz$
, and for (e,f)
$100$
μm thick water films vibrating at
$58\,\rm \rm kHz$
. The vibration amplitudes are
$a_0=$
5 μm (black symbols), 15 μm (red symbols) and 25 μm (blue symbols).
(a,c,e) Real and (b,d,f) imaginary parts of the rate factor
$\alpha$
for (a,b)
$100$
μm and (c,d)
$500$
μm thick ethanol films vibrating at
$20\,\rm \rm kHz$
, and for (e,f)
$100$
thick ethanol films vibrating at
$58\,\rm \rm kHz$
. The vibration amplitudes are
$a_0=$
5 μm (black symbols), 15 μm (red symbols) and 25 μm (blue symbols).
For a more global inspection of the unstable regimes of the liquid-film behaviour, the vibration amplitude
$a_0$
was varied beyond the three values represented in figures 4 and 5. The visualisation of the emerging results leads to three-dimensional surfaces of the disturbance growth rate
$\alpha _{\textit{real}}$
as a function of the disturbance wavenumber
$k$
and the vibration amplitude
$a_0$
. Figure 6 depicts two-dimensional projections of such surfaces into the (
$k,a_0$
) plane, revealing charts of the regimes of different stability of water and ethanol films of different thicknesses, vibrating at different frequencies. Depicted are only states with positive
$\alpha _{\textit{real}}$
, i.e. disturbances of the liquid-film surface growing in time. Unstable behaviour is found inside the regions on the nomograms bounded by the outermost, ‘tongue’-shaped lines. These outermost lines represent the neutral curve of the stability behaviour, i.e. states where
$\alpha _{\textit{real}} = 0$
. The colour information of the iso-lines of
$\alpha _{\textit{real}}$
for different values of
$\alpha _{\textit{real}}$
is additionally quantified by the values printed next to the lines. The charts additionally show the projections of the ‘ridges’ of the three-dimensional surfaces into the (
$k,a_0$
) plane as black dotted lines, marking the states of maximum film surface deformation growth rate for a given vibration amplitude
$a_0$
, i.e. the states which are most dangerous for the film in terms of its destabilisation. Furthermore, red solid lines mark the maximum growth rate
$\alpha _{\textit{real}}$
among all tongues for a given forcing amplitude
$a_0$
. The lines show that, for the smaller frequency of
$20\,\rm kHz$
, the global maxima coincide with the maxima in the leftmost tongue-shaped regime for these values of
$a_0$
. For the higher frequency of
$58\,\rm kHz$
, in contrast, with increasing amplitude
$a_0$
, the global maximum jumps from the leftmost to the neighbouring regime with higher wavenumber. The threshold value of
$a_0$
is approximately
$30\,\rm \unicode{x03BC} m$
for the case depicted.
Regimes of stable and unstable liquid-film deformation behaviour for varying non-dimensional wavenumber
$k$
and dimensional vibration amplitude
$a_0$
. The contour lines represent constant values of the dimensionless growth rate,
$\alpha _{\textit{real}}$
. Films with properties outside the tongue-shaped contours with
$\alpha _{\textit{real}} = 0$
are stable, inside they are unstable. (a,c,e) Water and (b,d,f) ethanol films. (a,b) Film thickness
$100\,\rm \unicode{x03BC} m$
, vibration frequency
$20\,\rm kHz$
; (c,d) film thickness
$500\,\rm \unicode{x03BC} m$
, vibration frequency
$20\,\rm kHz$
; (e,f) film thickness
$100\,\rm \unicode{x03BC} m$
, vibration frequency
$58\,\rm kHz$
. Black dashed lines mark the local maximum growth rate for a given
$a_0$
within each tongue. Red solid lines mark the global maximum growth rate for a given
$a_0$
across all tongues.
(a,d) Real and (b,e) imaginary parts and (c,f) complex plane of the rate factor
$\alpha$
for a
$500\,\rm \unicode{x03BC} m$
thick water film vibrating at
$20\,\rm kHz$
with an amplitude
$a_0$
of (a–c)
$5\,\rm \unicode{x03BC} m$
and (d–f)
$25\,\rm \unicode{x03BC} m$
, as of the present solutions and as per (5.1). The red dashed lines represent the solutions from (5.1).
To understand the imaginary component of the rate factor and the patterns formed in the
$\alpha _{\textit{imag}}$
versus wavenumber plots, we introduce figure 7. In this figure, we plot the real and imaginary parts of
$\alpha$
, together with the Gaussian planes of the rate factor, for
$500\,\rm \unicode{x03BC} m$
thick water films vibrating at
$20\,\rm kHz$
for two different vibration amplitudes. In the left-half Gaussian planes, where
$\alpha _{\textit{real}}\lt 0$
, the surface waves are recognised as stable and decay over time. We compare these solutions with the predictions by the dispersion relation for unforced, capillary surface waves in Rayleigh–Taylor instability (Chandrasekhar Reference Chandrasekhar1961). The dispersion relation quantifies the influences on the rate factor by the disturbance wavenumber, surface tension, liquid viscosity and density. Miles (Reference Miles1993) had already had the idea that Faraday waves could be understood as a vibrationally forced perturbation of capillary–gravity waves, having developed nonlinear theory for their pattern selection. Here, we present a solution that explicitly relates the dispersion relation of Faraday waves to the classical dispersion relation of capillary–gravity waves. For two layers of fluids, with a gas on top of a liquid, with the assumption that the modified wavenumber
$({\bar k}^2 + \bar {\alpha } / \nu )^{1/2}$
(overbars indicating dimensional quantities) is the same for both fluids and the liquid and gas densities are different by a factor
$O(10^3)$
, the dispersion relation assumes the form
\begin{equation} \alpha ^2 + 8 \frac {k^2 \nu }{\omega h^2} \alpha - 16 \frac {\nu ^2 k^4}{h^4 \omega ^2} \left ( -1 + \sqrt {1+\frac {h^2 \omega }{2 k^2 \nu } \alpha } \right ) + 4 \frac {gk}{h \omega ^2} + 4 \frac {k^3 \sigma }{h^3 \rho \omega ^2} = 0, \end{equation}
with non-dimensional
$\alpha$
and
$k$
. By solving this equation, the decay rate and oscillation frequency of unforced, viscous capillary surface waves are determined as the real and imaginary parts of the rate factor
$\alpha$
. Figure 7 depicts one solution to the Rayleigh–Taylor dispersion relation along with its complex conjugate, together with the solutions from the present dispersion relation, for two different forcing amplitudes
$a_0$
. The two solutions for the decay rate, i.e. in the case of stable film behaviour, agree perfectly at wavenumbers where no bifurcations are present, while the imaginary part of
$\alpha$
for the unforced, viscous capillary surface waves differs from our perturbed viscous surface wave solution in the stable regime for the larger
$a_0$
. The deviations decrease, however, with increasing wavenumber
$k$
, because the capillary and viscous forces dominate the inertial forces. This agrees with the findings by Rajchenbach & Clamond (Reference Rajchenbach and Clamond2015). The repeated pattern of
$\alpha _{\textit{imag}}$
emerges as a result of the complex conjugates and harmonics of solutions to the present dispersion relation.
In figure 7(a,d) we can begin examining the Hopf bifurcation. Starting at
$k=0$
, we see that the system is dynamically stable. As mentioned, the decay of the surface waves follows the solution of the stable Rayleigh–Taylor configuration. Moving to the right along the abscissa, we see the system go from aperiodically damped to oscillatory damped when the first Hopf bifurcation appears. Shortly following the bifurcation point, we see the system quickly turn into oscillatory unstable, as the top branch crosses into the positive plane. In figure 7(a), for the system described with a shaking amplitude of 5 μm, this occurs at a dimensionless wavenumber of
$k\approx 16$
. The system recovers to an aperiodically stable state as the wavenumber grows and again follows the solution of the stable Rayleigh–Taylor configuration between wavenumbers of 22 and 30. At
$k\approx 30$
we see the emergence of a second Hopf bifurcation towards oscillatory damped film behaviour. While this bifurcation selects a frequency, the film never becomes unstable. Figure 7(d) shows the same system previously discussed but for a larger shaking amplitude of 25 μm. Again, the system starts off as dynamically stable at
$k=0$
. The first two closed contours have grown when comparing with the 5 μm system. At this larger shaking amplitude, we see the second closed contour has reached the positive half-plane, and the film has become oscillatory unstable. Third and fourth bifurcations have also emerged; however, for the two associated ranges of
$\alpha _{\textit{real}}$
, the film surface deformation is oscillatory damped. Comparison of these solutions with the Rayleigh–Taylor dispersion relation shows the effects from viscosity, surface tension and gravity, which act to stabilise the surface against the input perturbations by shifting the ranges of
$\alpha _{\textit{real}}$
into the negative half-plane, where the surface waves are damped.
In figure 7(c,f), the vertical line drawn at
$\alpha _{\textit{real}}=0$
indicates the state of neutral stability, where the growth rate of the surface waves is zero. This corresponds to a neutral stability condition, implying neither damping nor amplification of the waves. We see that, at this threshold, the imaginary part of the rate factor can be non-zero. As the imaginary part of the rate factor corresponds to the response frequency of the surface waves, this means that a surface wave can perpetually exist while oscillating without growth or decay. In contrast,
$\alpha _{\textit{imag}}=0$
implies a non-oscillatory (monotonic) response of the surface deformation pattern.
Next, we examine the right-hand half-planes of figure 7(c,f) where
$\alpha _{\textit{real}} \gt 0$
, indicating the growth of surface waves. Within the unstable zone, integer values of
$\alpha _{\textit{imag}}$
are selected, indicating that only specific discrete frequencies emerge on an unstable film surface. These discrete values of frequencies align with the harmonic and subharmonic responses to the ultrasonic atomiser’s angular frequency. Expanding the sum in the formulation of the surface deformation in (3.9), we can deduce that the harmonic response of the film includes all even-valued integers of
$\alpha _{\textit{imag}}$
and the subharmonic response includes all odd-valued integers of
$\alpha _{\textit{imag}}$
. At the lowest perturbation amplitude,
$a_0=5$
μm, only the first tongue crosses the threshold for instability, shown in figures 4(a,b) and 5(a,b) for water and ethanol, respectively. This implies that only a subharmonic response is excited. In the unstable regime, the frequency is fixed for all wavenumbers that are contained within that unstable region. In contrast, at larger perturbation amplitude (
$a_0=25$
μm), both harmonic and subharmonic responses are excited. Although larger perturbation amplitudes raise the magnitude of the growth rate, they do not affect the frequency of the response, which is solely determined by the forcing frequency.
Depicting the regimes of unstable behaviour for films of a given liquid vibrating at a given frequency and amplitude, differing by their thickness, leads to the data in figure 8. The nomograms show that the neutral curves do not depend on the film thickness. Equation (4.6) for the quantity
$\eta _i$
confirms this. This behaviour is not unique to high frequencies and can also be observed analytically at frequencies as low as 100 Hz using the present method. Moreover, this independence of the neutral stability curve from film height has also been analytically confirmed across a range of scales between 0.1 and 10 mm. The disturbance growth rate, however, is proportional to
$h^{-1/2}$
, so that it increases with decreasing film thickness.
Regimes of stable and unstable water film deformation behaviour for varying dimensional wavenumber
$\bar k$
and dimensional vibration amplitude
$a_0$
, at a vibration frequency of
$40\,\rm kHz$
. The film thicknesses are (a) 100
$\rm \unicode{x03BC} m$
, (b) 300
$\rm \unicode{x03BC} m$
and (c) 500
$\rm \unicode{x03BC} m$
. (d) Neutral curves bounding the unstable regimes compiled together, showing the independence of the unstable regimes from the film thickness. The colour coding in (a–c) shows that the disturbance growth rate in the unstable regime is higher for thinner films.
Locations of states of operation of ultrasonic atomisers from Topp (Reference Topp1973) (black star), Sindayihebura (Reference Sindayihebura1995) (open square) and Ramisetty et al. (Reference Ramisetty, Pandit and Gogate2013) (black dotted line) in relation to the regimes of unstable water-film deformation behaviour as predicted by the present theory. Red solid lines mark the global maximum growth rate for a given
$a_0$
across all tongues.
5.2. Validation with other literature
In the literature, other authors reported studies of the presently investigated process. We first look at research investigating the Faraday instability of vibrating liquid films for spray formation, as reported in the visualisation studies by Topp (Reference Topp1973), Sindayihebura (Reference Sindayihebura1995) and Ramisetty et al. (Reference Ramisetty, Pandit and Gogate2013). In many papers on the subject, the sets of data about the liquid film and the state of operation of the atomiser are incomplete. Often the film thickness is unknown. From the top-view photographs of liquid vibrating films, however, the wavelength of the film surface deformation may be determined. The states of operation reported in the three papers referenced here are shown in figure 9, where the dimensional wavenumber on the abscissa
${\bar k} = 2 \pi / \lambda$
. Topp (Reference Topp1973) studied a water film vibrating at
$37\,\rm kHz$
, with the film thickness unknown. The wavelength of film surface deformation is found to be
$105\,\rm \unicode{x03BC} m$
, and the vibration amplitude
$a_0$
was
$4.5\,\rm \unicode{x03BC} m$
at the state of operation in Topp’s figure 2, where spray was formed. In his PhD thesis, Sindayihebura (Reference Sindayihebura1995) presented top-view photographs of vibrating water films on transducers of ultrasonic atomisers. The film thickness was between
$0.3$
and
$1\,\rm mm$
. For a vibration frequency of
$40.8\,\rm kHz$
, the film surface deformation wavelength is found to be
$144.92\,\rm \unicode{x03BC} m$
from his figure 5.9, with the vibration amplitude
$a_0$
of
$4.75\,\rm \unicode{x03BC} m$
. Ramisetty et al. (Reference Ramisetty, Pandit and Gogate2013) presented a top-view photograph of excellent quality of a vibrating water film on the transducer of an ultrasonic atomiser vibrating at
$20\,\rm kHz$
(see figure 2). The film surface deformation wavelength measured on the photograph is
$135.9\,\rm \unicode{x03BC} m$
, in contrast to the
$155\,\rm \unicode{x03BC} m$
reported by the authors. For a frequency of
$40\,\rm kHz$
, the wavelength reported is
$85\,\rm \unicode{x03BC} m$
. The film thickness is unknown. The vibration amplitude is unknown, but it was low enough to prevent spray formation, which would have disturbed taking the photograph. The findings from these three sources are depicted in figure 9 by symbols, where the nomogram axes are dimensional. We see that the theory reproduces the experimental findings in that the result by Topp (Reference Topp1973) lies in the unstable regime of the liquid film, the result by Sindayihebura (Reference Sindayihebura1995) lies in the stable regime and the state investigated by Ramisetty et al. (Reference Ramisetty, Pandit and Gogate2013) is represented by the black dotted line the figure, since the vibration amplitude is unknown. The results indicate that there exist values of the amplitude
$a_0$
leading to either stable or unstable film behaviour. For the stable state in figure 2, the amplitude had to be below
$7\,\rm \unicode{x03BC} m$
.
Kumar (Reference Kumar1996) studied the linear Faraday instability of viscous liquid films and presented regimes of stable and unstable film behaviour in
$(2 \pi / \lambda , a/g)$
maps. In figure 10 we compare the stability regimes presented in Kumar’s figure 1 for a glycerol–water mixture (
$\mu = 0.1246\,\rm Pa \, s$
,
$\rho = 1222\,\rm kg\,m^{-3}$
,
$\sigma = 67.6 \times 10^{-3}\,\rm N\,m^{-1}$
) with our theory. The film thickness is
$2\,\rm mm$
and the vibration frequency is
$60\,\rm Hz$
, i.e. not ultrasonic. The shapes of the neutral curves in the
$(2 \pi / \lambda , \textit{Fr})$
plane from the two works agree very well. This explicitly shows that we recover the solution of Kumar (Reference Kumar1996) and by extension that of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). In the range of small dimensional wavenumbers there are differences, which are due to the simplifications in the hyperbolic functions in our theory, which are accurate for sufficiently large wavenumber
$k$
only.
Regimes of unstable glycerol–water film deformation behaviour for a square-wave pattern with varying Froude number and dimensional disturbance wavenumber, compared with Kumar (Reference Kumar1996). The regimes inside the tongue-shaped boundaries are characterised by unstable film deformation behaviour. The red lines mark Kumar’s results, the coloured regimes are the present results. In the shaded region, our solution to the dispersion relation is less accurate due to the approximated form of
$\eta _i$
.
Sindayihebura (Reference Sindayihebura1995) and Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) analysed theoretically the Faraday instability of a viscous liquid film, assuming a surface-wave pattern governed by functions of the form
$\cos n \theta \,J_n(kr)$
(polar angle
$\theta$
in the horizontal plane, Bessel function
$J_n$
), rather than the square-shaped pattern found experimentally in resonance. Figure 11 compares the stability regimes of water films at constant Weber and Reynolds numbers for four different values of the non-dimensional vibration frequency
$\varOmega$
. The figure displays the unstable regimes at the lowest wavenumbers for each of the four different frequencies, both as reported by Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998) in their figure 5 and as obtained by the present analysis. It is seen that the effects from the different approaches of the two works and from the errors are not very large. The principal shapes of the regimes of unstable liquid-film behaviour are very similar despite Sindayihebura and Bolle’s representation of the square planform with cosine and Bessel functions. In principle, a Bessel–Fourier expansion may be used to approximate the square planform which was observed experimentally; however, this is not the most efficient way of representing the observed planform, and is better suited for patterns with point symmetry.
Regimes of unstable water-film deformation behaviour for varying non-dimensional vibration frequency at
$\textit{We} = 0.13$
and
$\textit{Re} = 66$
as of the present solutions, compared with Sindayihebura & Bolle (Reference Sindayihebura and Bolle1998). The regimes with the minimal Froude numbers increasing with the corresponding wavenumber correspond to
$\varOmega =$
$1268$
,
$2538$
,
$3172$
and
$4123$
.
5.3. Minimal required vibration amplitude and global maximum of the growth rate
For forming a spray by destabilisation of a liquid film on the transducer surface of an ultrasonic atomiser, for a given liquid, film thickness and vibration frequency, a minimum vibration amplitude must be exceeded. Those minimum values of
$a_0$
required are given by the tip points on the bottoms of the tongue-shaped regimes in the instability charts, such as figures 6 and 8. A non-dimensional correlation for those values covering all the cases presently studied was developed, using dimensional analysis. To determine the coefficients for the dimensional analysis, Matlab’s built-in optimisation function ‘fminsearch’ was used to minimise the least squares error in the linear regression. The resulting coefficients were then rounded to simple fractions, which are the values presented. The dimensional parameters entering the analysis, together with the critical vibration amplitude
$a_{0,cr}$
as the target quantity, the film thickness
$h$
and the gravitational acceleration
$g$
, are the liquid density
$\rho$
, dynamic viscosity
$\mu$
and surface tension against air
$\sigma$
, and the vibration frequency
$f = \omega / 2 \pi$
. The ranges of the latter dimensional parameters covered by the analysis are
$789\,\rm kg\,m^{-3} \leqslant \rho \leqslant 1224\,\rm kg\,m^{-3}$
,
$0.59\,\rm mPa\, s \leqslant \mu \leqslant 130\,\rm mPa \, s$
and
$22 \times 10^{-3}\,\rm N\,m^{-1} \leqslant \sigma \leqslant 72 \times 10^{-3}\,\rm N\,m^{-1}$
. The vibration frequency varied in the range
$60\,{\rm Hz} \leqslant f \leqslant 1\,{\rm MHz}$
. With three dimensions involved, the result is a correlation of the non-dimensional vibration amplitude
$a_{0,cr} / h$
with the three quantities
$\varOmega$
,
$\textit{We}$
and
$\textit{Re}$
. The correlation reading
\begin{equation} \frac {a_{0,cr,SH}}{h} = \frac {4}{3} \varOmega ^{-2/5}\, \textit{We}^{1/5}\, \textit{Re}^{-4/5} \end{equation}
represents the tip of the subharmonic, leftmost tongue in the stability charts and includes the stabilising effects of the surface tension and of the dynamic viscosity of the liquid. The correlation can be rewritten in the form
\begin{equation} \frac {a_{0,\textit{cr},\textit{SH}}}{h} = \frac {4}{3} \left ( \frac {Oh^2}{\omega ^*} \right )^{2/5}\!, \end{equation}
where
$Oh = \mu / (\sigma h \rho )^{1/2}$
and
$\omega ^* = \omega (\rho h^3 / \sigma )^{1/2}$
. Inspection of the correlation (5.2) shows that the gravitational acceleration cancels, and that, in a dimensional representation of
$a_{0,cr}$
, the film thickness also cancels. The film thickness
$h$
in (5.3), therefore, plays the role of a scaling factor only. These results represent the smallest vibration amplitude required to destabilise the liquid film for the purpose of spray formation. Only data points with negligible deviation of the approximated dispersion relation from the exact form are included. The data are shown in figure 12(a), together with Kumar’s results which align well with the present theory. At large values of
$a_{0,cr} / h$
, however, the data deviate, which is due to the small wavenumbers associated with the states of operation shown. This is seen in figure 10, where, at low
$k$
, the present correlation underpredicts that by Kumar (Reference Kumar1996).
Universal representation of the minimal vibration amplitude needed for unstable liquid-film behaviour (a) on the leftmost (subharmonic) tongue and (b) on the second tongue from the left (harmonic). (c) Non-dimensional vibration amplitude at the transition of the maximum growth rate from the first to the second tongue.
An analogous analysis of the tip of the second, harmonic tongue in the stability charts yields the correlation
\begin{equation} \frac {a_{0,cr,H}}{h} = \frac {25}{4} \left ( \frac {Oh^{13}}{{\omega ^*}^9} \right )^{1/20}\!. \end{equation}
With increasing forcing amplitude and frequency, furthermore, the global maximum of the disturbance growth rate, which is located in the subharmonic regime, i.e. in the leftmost tongue of the stability chart, for small frequencies, jumps to the harmonic regime to the right of the former regime. The critical forcing amplitude characterising this transition is denoted
$a_{0,cr,t}$
. For this quantity, the correlation
\begin{equation} \frac {a_{0,cr,t}}{h} = \frac {77}{10} \left ( \frac {Oh}{{\omega ^*}^3} \right )^{1/5} \end{equation}
was found. This is an important relation for the operation of ultrasonic atomisers in view of the dominant wavelength of the liquid-film surface wave patterns, which influences the drop size formed by the destabilisation of the film by the vibrations. The scaling factor
$h$
drops out from (5.4) and (5.5) also.
6. Conclusions
A linear analysis of the Faraday instability for a liquid film on a vibrating substrate is carried out with the aim of studying the influences of the vibration frequency, vibration amplitude, liquid-film thickness and liquid material parameters on the unstable regimes of the wavenumber, and growth or damping rates of the square-patterned surface deformations. The present method is applicable across a broad frequency range, from low frequencies (
${\lt} 100\, \rm Hz$
) to high frequencies (tested up to 1 MHz). Limitations in the accuracy of the dispersion relation at low disturbance wavenumbers are only due to a simplification of hyperbolic functions in a term. The results quantify not only the ranges of vibration wavenumbers and amplitudes needed to destabilise the liquid film, i.e. suitable for ultrasonic atomisation of the film for the purpose of spray formation, but also the growth rate in time and the related vibration frequency inside the unstable regimes of the film-surface deformation. The frequency of the oscillations is shown to relate to the viscous dispersion of unforced capillary waves. We examine the Hopf bifurcations, and observe how the system transitions from aperiodically damped to oscillatory damped, and finally we see ranges where the film becomes oscillatory unstable. In this analysis, we can explicitly see how the surface tension, viscosity and gravity act as restoring forces to stabilise the perturbed liquid film. For a given liquid, higher vibration frequencies produce unstable film surface waves with shorter wavelengths. Increasing vibration amplitudes and decreasing film thickness produce higher growth rates of the surface deformation in time. Increasing surface tension of the liquid against the ambient air reduces the deformation growth rate. With increasing vibration amplitude
$a_0$
, the global maximum of the deformation growth rate occurs in different unstable regimes in the
$(a_0,k)$
plane. In addition, it is shown that a jump of the global maximum of the deformation growth rate from the first, subharmonic tongue to the neighbouring, harmonic tongue occurs at different values of
$a_0/h$
. The value of
$a_0/h$
where the tongue transition occurs depends on the fluid properties and the operating parameters of the vibrating surface. Non-dimensional relations are found by dimensional analysis to represent the minimum vibration amplitudes required to destabilise the film in the subharmonic and harmonic regimes, and the vibration amplitude for the global maximum disturbance growth rate transition from the subharmonic to the harmonic regime. The relations show the stabilising roles of the liquid density and dynamic viscosity, as well as the surface tension against ambient air.
Acknowledgements
Financial support of this research project by the German Research Foundation (DFG) through project number BR 1046/15-1 in the framework of the Special Research Program SPP 2419 is gratefully acknowledged. The authors gratefully acknowledge the support by Dr D. Zrnić in solving the dispersion relation.
Declaration of interests
The authors report no conflicts of interest.
Open access
