2.1 Governing equations

Fluid 1 is hydrodynamically active and satisfies the continuity and Navier–Stokes equations. Moreover, as it is a perfect conductor, the electrostatic potential everywhere in the domain equals the bottom wall potential. These equations are given by

Here, $\unicode[STIX]{x1D64F}_{1}=-p_{1}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}_{1}(\unicode[STIX]{x1D735}\boldsymbol{v}_{1}+\unicode[STIX]{x1D735}\boldsymbol{v}_{1}^{\text{T}})$ , $\boldsymbol{i}_{z}$ is the unit vector along the positive $z$ direction and $\unicode[STIX]{x1D644}$ is the identity tensor. For fluid 2, the electric field, $\boldsymbol{E}_{2}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{2}$ , satisfies Gauss’s law,

The interface speed ( ${\mathcal{U}}$ ), the unit normal vector ( $\boldsymbol{n}$ ) and the unit tangent vector ( $\boldsymbol{t}$ ) are given by

At the fluid–fluid interface, $\boldsymbol{v}_{1}\boldsymbol{\cdot }\boldsymbol{n}-{\mathcal{U}}=0$ must hold since it is a material surface, while the interfacial stress balances are given by

where $\unicode[STIX]{x1D648}_{2}=\unicode[STIX]{x1D700}_{0}\unicode[STIX]{x1D700}_{2}[\boldsymbol{E}_{2}\boldsymbol{E}_{2}-1/2(\boldsymbol{E}_{2}\boldsymbol{\cdot }\boldsymbol{E}_{2})\unicode[STIX]{x1D644}]$ is the Maxwell stress exerted by fluid 2 (Saville Reference Saville1997). We now derive the long-wave model based on the separation of length scales.

2.2 Long-wave model

The governing equations are non-dimensionalized using the following scales, denoted by the subscript ‘ $c$ ’:

Here, $\unicode[STIX]{x1D716}$ is a small parameter defined as $\unicode[STIX]{x1D716}=H/\unicode[STIX]{x1D6EC}$ , known as the ‘film parameter’. The characteristic velocity scale, $U$ , is determined from the forcing frequency and given by $U=\unicode[STIX]{x1D714}H/2\unicode[STIX]{x03C0}$ or $U=Hf$ . Using the above scales, the non-dimensional flow equations, retaining terms up to $O(\unicode[STIX]{x1D716})$ and dropping higher-order terms, are obtained as

and

To obtain the above governing equations, $Re$ is taken to be at most $O(1)$ , while $G$ to be at least $O(1)$ . The no-slip and no-penetration conditions are given as $u_{1}(0)=w_{1}(0)=0$ . Now, turning to the electrostatic equations, we obtain $\unicode[STIX]{x1D6F7}_{1}=\cos (2\unicode[STIX]{x03C0}t)$ . For fluid 2, we obtain $\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F7}_{2}/\unicode[STIX]{x2202}z^{2}=0$ from (2.1b ), subject to the boundary conditions, $\unicode[STIX]{x1D6F7}_{2}(h)=\cos (2\unicode[STIX]{x03C0}t)$ and $\unicode[STIX]{x1D6F7}_{2}(\unicode[STIX]{x1D6FD})=0$ , where $\unicode[STIX]{x1D6FD}=H_{w}/H$ . Solving for $\unicode[STIX]{x1D6F7}_{2}$ , we obtain $\unicode[STIX]{x1D6F7}_{2}=(\unicode[STIX]{x1D6FD}-z)/(\unicode[STIX]{x1D6FD}-h)\cos (2\unicode[STIX]{x03C0}t)$ .

The tangential and normal stress balance at the interface ( $z=h$ ) gives

In the above, $Ca$ is taken to be at most $O(\unicode[STIX]{x1D716}^{2})$ , which corresponds to the large-surface-tension limit, true for fluids such as water or liquid metal conductors. In addition, $E$ is taken to be at least $O(1)$ . Using a rough order of magnitude analysis for the water–air system, it can be verified that all of the approximations made in the model are practical. Taking $\unicode[STIX]{x1D70C}_{1}\sim O(10^{3})$ , $\unicode[STIX]{x1D707}_{1}\sim O(10^{-3})$ , $\unicode[STIX]{x1D6FE}\sim O(10^{-1})$ , $f\sim O(1)$ , $H\sim O(10^{-3})$ , $\unicode[STIX]{x1D700}_{2}\sim O(1)$ , $\unicode[STIX]{x1D700}_{0}\sim O(10^{-11})$ and $A\sim O(10^{3})$ , it can be verified that all parameters ( $Re,G,Ca$ and $E$ ) satisfy the order of magnitude assignments made in the model.

By integrating the continuity equation and using the kinematic condition, we obtain $\unicode[STIX]{x2202}h/\unicode[STIX]{x2202}t=-\unicode[STIX]{x2202}q/\unicode[STIX]{x2202}x$ , where $q$ is the flow rate given by $q=\int _{0}^{h}u_{1}\,\text{d}z$ . By integrating the vertical component of momentum balance and eliminating pressure using (2.6) and (2.7), we obtain

We now use the weighted residual technique as described by Ruyer-Quil & Manneville (Reference Ruyer-Quil and Manneville2002) to derive the set of reduced-order evolution equations that govern the interface position, $h(x,t)$ , and flow rate, $q(x,t)$ . Towards this end, we first decompose the horizontal velocity component as follows:

Here, $\hat{u} _{1}$ is the $O(1)$ contribution to the velocity while $\tilde{u} _{1}$ denotes the $O(\unicode[STIX]{x1D716})$ correction. The vertical component of velocity, $w_{1}$ , as shown in (2.9b ), is then obtained by integrating the continuity equation. The leading-order velocity, $\hat{u} _{1}$ , is chosen to be locally parabolic at each horizontal position. This assumption is justified for moderate Reynolds numbers. Thus, $\hat{u} _{1}$ is determined using the following equations:

Here, $K_{u}$ is introduced so that the leading-order velocity profile is locally parabolic, and is obtained in terms of the flow rate, $q$ , using the integral constraint in (2.10d ). Next, we substitute for velocity as given by (2.9) in (2.8) and neglect all terms of $O(\unicode[STIX]{x1D716}\tilde{u} _{1})=O(\unicode[STIX]{x1D716}^{2})$ or smaller. Thus, the velocity corrections that appear in the inertial terms (left-hand side of (2.8)) are ignored. The only term that still contains the $O(\unicode[STIX]{x1D716})$ term, $\tilde{u} _{1}$ , is due to viscous diffusion. This term is eliminated by taking a weighted integral of (2.8) with a suitable weight function. We use the Galerkin method, wherein the weight function has the same functional form as the leading-order velocity, as it is the most efficient choice of weight function (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2002). The weight function, $F$ , is then defined as follows:

Using the weighted residual strategy, the diffusion term now gives $\int _{0}^{h}F(\unicode[STIX]{x2202}^{2}u_{1}/\unicode[STIX]{x2202}z^{2})\,\text{d}z=q$ upon integration by parts. The WRIBL equation obtained from (2.8) is then given by

Equation (2.12) represents the long-wave WRIBL model consistent at $O(\unicode[STIX]{x1D716})$ . Substituting for velocity and $F$ obtained from (2.10) and (2.11), we obtain the evolution equations in terms of $h$ and $q$ alone. Once the reduced-order model is obtained, we renormalize the system by setting $\unicode[STIX]{x1D716}=1$ , as done by Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2013). This gives the final set of evolution equations as

In (2.13b ), the left-hand side consists of inertial terms. These inertial contributions were also obtained using an alternative method of polynomial expansion by Rojas et al. (Reference Rojas, Argentina, Erda and Tirapegui2010) and also by Bestehorn (Reference Bestehorn2013) for the mechanical Faraday problem. It should be noted that Rojas et al. (Reference Rojas, Argentina, Erda and Tirapegui2010) derived an $\unicode[STIX]{x1D716}^{2}$ model and therefore retained additional effects of lateral diffusion as well. On the right-hand side of (2.13b ), the first term is due to transverse viscous diffusion, while the terms inside the bracket correspond to electrical, surface tension and gravitational effects respectively. It should be noted that when $E=0$ , the model reduces to an inertial lubrication model for inertial surface-gravity waves on thin liquid films. Under this limit, the presence of the gas-bounding wall is not felt by the fluid, as fluid 2 is assumed to be hydrodynamically passive.

The linear stability analysis of (2.13) is carried out using the method of modal decomposition. In the base state, the fluid is assumed to be stationary with a flat interface. The variables are perturbed around the base state as $h=1+h^{\star }\text{e}^{\text{i}kx}$ and $q=q^{\star }\text{e}^{\text{i}kx}$ , where the starred variables denote the perturbation amplitudes and $k$ is the spatial wavenumber. This leads to the following damped Mathieu equation:

In the above, $\unicode[STIX]{x1D714}_{0}$ is the natural frequency of the interface. It may be noted that the natural frequency of a thin film is modified in the presence of an electrostatic forcing (cf. the additional term due to the Maxwell stress). It can be observed that in (2.14), for a gravitationally unstable system, $G<0$ , and that the electric term is of the same sign as the gravity term. This observation forecasts the result that periodic electrostatic forcing cannot stabilize an RT unstable configuration of a conducting liquid, such as water, lying above a dielectric gas, such as air. Moreover, the Mathieu equation has a forcing frequency $4\unicode[STIX]{x03C0}$ compared with the imposed potential frequency $2\unicode[STIX]{x03C0}$ . This is due to the quadratic dependence of the Maxwell stress on the electric field. A stability analysis of the Mathieu equation (2.14) is carried out using standard Floquet theory (Nayfeh Reference Nayfeh1981).

For nonlinear calculations, an initial wavy disturbance is imposed on the interface position and equations (2.13) are solved on a computational domain corresponding to one wavelength of the imposed perturbation. The initial conditions are $h(x,0)=1+0.05\,\cos (kx)$ and $q(x,0)=0$ unless otherwise specified. We use the Fourier spectral collocation technique to obtain numerical solutions (Labrosse Reference Labrosse2011). It ensures high-order spatial resolution and that periodic boundary conditions are satisfied. The number of grid points, $N$ , was chosen to be 120 for all calculations, and varied between 200 and 300 for sliding-related calculations. The time integration is performed using the adaptive time-stepping NDSolve subroutine in Mathematica^® v.10.3.

3.1 Faraday instability in a gravitationally stable configuration

We begin with the gravitationally stable configuration. In figure 2(a), the neutral stability curve for such a case is plotted for a thin film of 0.25 mm. The dashed line denotes the neutral curve obtained in the absence of inertia ( $Re=0$ ). It is seen that the system becomes unstable beyond a critical threshold potential to a long-wave instability. This corresponds to the ‘pillaring’ instability reported by Wu & Chou (Reference Wu and Chou2003) and Roberts & Kumar (Reference Roberts and Kumar2009). In the presence of fluid inertia (solid), additional stability tongues corresponding to alternating subharmonic and harmonic Faraday responses are observed at higher threshold potentials and larger wavenumbers. In figure 2(b), the neutral stability curve is plotted for a thicker film of 0.75 mm. For this case, the instability threshold for the Faraday modes is lower than that for the pillaring mode. Thus, Faraday modes are expected to be observed in an experiment at the onset of instability for thicker films, and more so as the liquid thickness increases, i.e. as the fluid becomes more inertial. In this work, we restrict our discussion to such inertial thin films of $O(1~\text{mm})$ and proceed to discuss their interfacial dynamics.

Figure 2. (a) The neutral stability curve (solid) shows the long-wave ‘pillaring’ instability along with alternating subharmonic (SH) and harmonic (H) tongues characteristic of Faraday resonance. In the absence of inertia (dashed), only the ‘pillaring’ mode exists; $H=0.25~\text{mm}$ , $\unicode[STIX]{x1D6FD}=5$ , $f=3~\text{Hz}$ . (b) The neutral stability curve with $H=0.75~\text{mm}$ , $\unicode[STIX]{x1D6FD}=5$ , $f=3~\text{Hz}$ ; the markers denote the parameters chosen for the nonlinear simulations.

Figure 3. (a) The linear (dashed) and nonlinear (solid) temporal evolution of $h_{max}$ for the pillaring mode at ▴ in figure 2(b) ( $k=0.1$ , $A=8.7~\text{kV}$ ). (b) The temporal evolution of $h_{min}$ for lower liquid holdup $(k=0.1,A=28.5~\text{kV},\unicode[STIX]{x1D6FD}=10)$ ; see movie 1 available at https://doi.org/10.1017/jfm.2018.682 for the interface evolution.

It is known that the solution to the Mathieu equation can become unbounded in two ways, namely (i) a non-monotonic non-periodic unbounded solution and (ii) a periodic solution with an exponentially increasing amplitude (Nayfeh Reference Nayfeh1981). We show that these correspond to the pillaring-type and Faraday modes respectively. In figure 3(a), the temporal evolution of the maximum position of the interface, $h_{max}$ , corresponding to the pillaring mode ( $\blacktriangle$ in figure 2 b) is plotted. The dashed curve denotes the linear Mathieu response while the solid curve denotes the nonlinear evolution. It can be seen that the Mathieu response is of the first kind and the corresponding nonlinear evolution of the interface shows a subcritical behaviour, with the interface approaching the top wall. As inertial films are typically of $O(1~\text{mm})$ , we do not consider long-range molecular forces, which become important only when the film thins to very small length scales of $O(100~\text{nm})$ (Israelachvili Reference Israelachvili2011). These forces are, however, of significance if one is interested in studying the spinodal dewetting and contact line dynamics. In the presence of these forces, pillar formation is arrested near the wall and ‘quasi-steady-state oscillations’ are observed, as shown by Roberts & Kumar (Reference Roberts and Kumar2009). In figure 3(b), a lower liquid–gas holdup ( $\unicode[STIX]{x1D6FD}=10$ ) is considered and the minimum position of the interface, $h_{min}$ , as a function of time is plotted. Here, the system again exhibits a subcritical behaviour, with the interface approaching the bottom wall (also see movie 1). The nonlinear emergence of the higher unstable Faraday modes, along with the viscosity of the film, considerably slows down the approach of the interface towards the bottom wall.

Figure 4. (a) The linear (dashed) and nonlinear (solid) temporal evolution of the interface position at $x=0$ , $h_{0}$ for the parameters corresponding to the markers in figure 2(b). (a) The subharmonic response at ●  $(k=0.19,5.6~\text{kV})$ . (b) The harmonic response at ▪  $(k=0.33,8.1~\text{kV})$ . Also see movie 2 for the interface evolution.

Figure 5. The effect of the forcing amplitude on the subharmonic response. (a) Large liquid–gas holdup: dashed $(A=6.7~\text{kV})$ , solid $(A=7.6~\text{kV})$ ; $k=0.22,\unicode[STIX]{x1D6FD}=5$ . (b) Low liquid–gas holdup: solid $(A=22~\text{kV})$ , dotted $(A=26~\text{kV})$ ; $k=0.22,\unicode[STIX]{x1D6FD}=10$ .

In figure 4(a,b), the linear (dashed) and nonlinear (solid) responses corresponding to subharmonic and harmonic Faraday modes are plotted (also see movie 2). It can be seen that the linear Mathieu response is oscillatory with an exponentially increasing amplitude, corresponding to the second kind. From figure 4(a), it is evident that the subharmonic response of the interface saturates nonlinearly to standing waves, isosynchronous with the forcing frequency. This is because the forcing ‘felt’ by the interface due to the electrostatic Maxwell stress has a quadratic dependence on the applied potential. Similarly, the harmonic response has a frequency that is twice that of the forcing frequency, as is evident from figure 4(b), and the interface oscillates twice every forcing cycle. The amplitude of the standing waves increases with an increase in the forcing amplitude, and the interface may approach either wall depending on the liquid–gas holdup. Figure 5(a) shows that for a large liquid–gas holdup ( $\unicode[STIX]{x1D6FD}=5$ ), the wave amplitude increases with an increase in the applied potential, and the interface spikes up towards the top wall. The inset of figure 5(b) is a plot of $(h_{max}-1)$ of the standing wave solutions as a function of $A^{2}-A_{c}^{2}$ . Here, $A^{2}-A_{c}^{2}$ is the bifurcation parameter, while $(h_{max}-1)$ acts as an ‘order parameter’ of the standing wave solutions (Cross & Greenside Reference Cross and Greenside2009). Simulations are carried out with various initial conditions of the form $h(x,0)=1+\unicode[STIX]{x1D712}\,\cos (kx)$ and $q(x,0)=0$ , wherein $\unicode[STIX]{x1D712}$ is varied from 0.005 to 0.9 and therefore also includes finite-amplitude perturbations. It is clear that the interface saturates to standing waves of amplitude ${\sim}0.9$ (cf. branch 2 in the inset) when $A^{2}>A_{c}^{2}$ , and that it is a subcritical bifurcation exhibiting a double hysteresis. For $A^{2}<A_{c}^{2}$ , very small disturbances decay to zero (cf. diamonds), consistent with the linear theory, while finite-amplitude disturbances may saturate to either branch 1 (cf. circles) or branch 2 (cf. squares) depending on the initial amplitude of perturbation, $\unicode[STIX]{x1D712}$ . Dynamic numeric calculations can saturate only to stable solutions. The unstable solutions are expected to emerge backward from the origin connecting the two stable branches (1 and 2). For lower liquid–gas holdup $(\unicode[STIX]{x1D6FD}=10)$ , the interface approaches the bottom wall with an increase in the forcing amplitude (cf. figure 5 b).

Figure 6. (a) Sliding of the interface as it approaches the wall; $t=0$ (dotted), $t=0.3\times 10^{5}$ (dashed), $t=2\times 10^{5}$ (dot-dashed), $t=5.8\times 10^{5}$ (solid). (b) The temporal evolution of $h_{min}$ . The markers denote the times corresponding to the interface profiles plotted in (a); $\unicode[STIX]{x1D6FD}=10,\unicode[STIX]{x1D714}=0,A=26~\text{kV},k=0.2.$ Also see movie 3 for the interface evolution.

From (2.13), we observe that in the limit of $\unicode[STIX]{x1D6FD}\gg h$ and low frequency ( $\unicode[STIX]{x1D714}\rightarrow 0$ ), the system reduces to an RT problem with an effective gravity, $\bar{G}$ , given by $\bar{G}=2E/\unicode[STIX]{x1D6FD}^{3}-G$ . In figure 6(a), we plot the interface profile with time for such a case, where $\unicode[STIX]{x1D6FD}=10$ and $\unicode[STIX]{x1D714}=0$ . The corresponding evolution of $h_{min}$ is plotted in figure 6(b). The markers denote the times corresponding to the interface profiles in figure 6(a). Since $\unicode[STIX]{x1D714}=0$ , the time on the abscissa in figure 6(b) is scaled with the capillary time, $\unicode[STIX]{x1D707}H/\unicode[STIX]{x1D6FE}$ . It is seen in figure 6(a) that the primary trough of the interface undergoes a buckling to form secondary troughs as it approaches the wall. This buckled quasi-steady configuration then continues to evolve at a very slow rate and subsequently exhibits ‘sliding’ (also see movie 3) as a consequence of a symmetry-breaking instability, wherein the interface translates parallel to the wall (cf. Lister, Morrison & Rallison Reference Lister, Morrison and Rallison2006; Glasner Reference Glasner2007; Lister, Rallison & Rees Reference Lister, Rallison and Rees2010; Dietze & Ruyer-Quil Reference Dietze and Ruyer-Quil2015). The necessary asymmetric perturbation that triggers the sliding motion is provided by numerical truncation errors, as detailed by Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015). In an experimental setting, it is impossible to shield the quasi-steady interface from random infinitesimal perturbations, which include symmetric as well as asymmetric components. The presence of any such asymmetric contribution will result in the interface exhibiting sliding dynamics. The onset of sliding coincides with the sharp decrease of $h_{min}$ with time at around $t=5\times 10^{5}$ in figure 6(b).

3.2 Unstable gravitational configuration

It is known that mechanical Faraday forcing as well as parallel electrical forcing is capable of stabilizing an RT unstable interface (Cimpeanu et al. Reference Cimpeanu, Papageorgiou and Petropoulos2014; Sterman-Cohen et al. Reference Sterman-Cohen, Bestehorn and Oron2017). Here, we explore the possibility of stabilizing an RT unstable interface using periodic perpendicular electrostatic forcing. In figure 7(a), the neutral stability curve for such a case is plotted (solid curve). The dashed vertical line corresponds to the RT cutoff wavenumber in the absence of forcing. Electrostatic forcing only renders the system more unstable, shifting the cutoff wavenumber to the right and increasing the range of RT unstable wavenumbers. In addition, it also results in the finite-wavelength subharmonic and harmonic Faraday instability for large-amplitude forcings. These Faraday modes behave very similarly to those discussed in the earlier section on the gravitationally stable configuration. Thus, as forecast earlier in § 2.1, the thin-film model shows that an RT unstable configuration of a heavy perfect conducting liquid over a lighter perfect dielectric cannot be stabilized by periodic electrostatic forcing. In figure 7(b), $h_{min}$ is plotted with time for the RT unstable configuration with forcing frequencies of 0.2 Hz (solid) and 0.04 Hz (dot-dashed), while the dashed curve corresponds to pure RT instability in the absence of forcing. In order to make a meaningful comparison, the time on the abscissa in figure 7(b) is again scaled with the capillary time, $\unicode[STIX]{x1D707}H/\unicode[STIX]{x1D6FE}$ . It is seen that around $t\sim 30\times 10^{5}$ , all $h_{min}$ curves exhibit a sudden change in their behaviour. This corresponds to the onset of sliding, and in the presence of periodic electrostatic forcing, the interface translates parallel to the wall, exhibiting an ‘earthworm-like’ motion (see movie 4).

Figure 7. (a) The neutral stability curve for an RT unstable interface with AC forcing $(\unicode[STIX]{x1D6FD}=5,f=3~\text{Hz})$ . The dashed vertical line is the RT cutoff wavenumber in the absence of forcing. (b) The temporal evolution of $h_{min}$ : dashed, pure RT ( $k=0.2,A=0$ ); solid, RT with forcing $(f=0.2~\text{Hz},A=2~\text{kV},k=0.2)$ ; dot-dashed, RT with forcing $(f=0.05~\text{Hz},A=2~\text{kV},k=0.2)$ . Also see movie 4.