Interface coupling effect and multi-mode Faraday instabilities in a three-layer fluid system

Abstract

We investigate the Faraday instabilities of a three-layer fluid system in a cylindrical container containing low-viscosity liquid metal, sodium hydroxide solution and air by establishing the Mathieu equations with considering the viscous model derived by Labrador et al. (J. Phys.: Conf. Ser., vol. 2090, 2021, 012088). The Floquet analysis, asymptotic analysis, direct numerical simulation and experimental method are adopted in the present study. We obtain the dispersion relations and critical oscillation amplitudes of zigzag and varicose modes from the analysis of the Mathieu equations, which agree well with the experimental result. Furthermore, considering the coupling strength of two interfaces, besides zigzag and varicose modes, we find a beating instability mode that contains two primary frequencies, with its average frequency equalling half of the external excitation frequency in the strongly coupled system. In the weakly coupled system, the $A$-interface instability, $B$-interface instability and $A$&$B$-interface instability are defined. Finally, we obtain a critical wavenumber $k_c$ that can determine the transition from zigzag or varicose modes to the corresponding $A$-interface or $B$-interface instability.

1. Introduction

The Faraday wave is a well-known phenomenon in fluid mechanics that has received considerable attention (Benjamin & Ursell 1954; Douady 1990; Milner 1991; Kumar 1996; Miles 1999; Rajchenbach, Leroux & Clamond 2011; Shao et al. 2021c). It is a nonlinear standing wave pattern generated by parametric excitation at the fluid interface (Faraday 1831). The primary frequency of these waves is equal to half the excitation frequency. Applications of the Faraday wave can be found in industrial areas such as self-assembly of particles (Chen et al. 2014) and cells (Guex et al. 2021), controlling chemical reactions (Hwang et al. 2020), and so on. To date, studies have been conducted on container shapes (Henderson & Miles 1991; Milner 1991; Umeki 1991; Rajchenbach et al. 2011; Batson, Zoueshtiagh & Narayanan 2013, 2015; Laroche et al. 2019), multiple physical fields (Paul & Kumar 2007; Zhao, Tang & Liu 2018; Ward, Matsumoto & Narayanan 2019a; Brosius et al. 2021), contact lines (Shao et al. 2021c; Wilson et al. 2022) and chaos (Guowei & Jiachun 1997). Corresponding phenomena, such as stripes, squares, hexagons and pentagonal stars, have been observed (Cross & Hohenberg 1993; Rajchenbach, Clamond & Leroux 2013). However, most of the previous studies have focused on the classical Faraday wave system with a single interface, such as a gas–liquid interface or a liquid–liquid interface. There is a significant lack of studies on Faraday waves with multiple interfaces.

More recently, interest in the two-interface Faraday instability has been revived by the work of Pucci et al. (2011) and Pucci, Amar & Couder (2013), which has shown the possibility of interesting phenomena due to the coupling between the two interfaces. By introducing a flexible boundary effect, Pucci et al. (2011, 2013) extended the classical Faraday instability to scenarios with two deformable interfaces (liquid–liquid and liquid–gas). The response of a droplet, suspended in a cavity filled with a viscous liquid of limited depth, to an external vibration has been investigated. It is found that an equilibrium has been reached between the radiation pressure exerted by Faraday waves on the borders and their capillary response. However, the coupling effects between interfaces and the influence of viscosity were ignored when estimating the dispersion relation. Pototsky & Bestehorn (2016), Pototsky et al. (2004), Pototsky, Oron & Bestehorn (2019) and Bestehorn & Pototsky (2016) performed a detailed numerical simulation of Faraday waves in a gas–liquid–liquid viscous system with strongly coupled interfaces. Zigzag and varicose modes of Faraday instability in a three-layer fluid system were identified from both experimental results and Floquet analysis using complete hydrodynamic equations. Regarding the zigzag mode, also called the barotropic mode, interface $A$ separating fluid 1 from fluid 2 oscillates in phase with surface $B$ separating fluid 2 from fluid 3, as shown in figure 1(b). Similarly, as shown in figure 1(c), the varicose mode involves oscillations in which the two interfaces exhibit out-of-phase behaviour. However, there is still a need for further clarification of critical information regarding the correlation between instability modes and the dispersion relation.

Figure 1. (a) Schematic diagram of the geometric model. (b) Zigzag mode. (c) Varicose mode.

From a theoretical aspect, Benjamin & Ursell (1954) first adopted linear stability analysis as an essential tool for studying Faraday waves to obtain the Mathieu equation for inviscid fluids. The critical conditions for Faraday instability were found on the boundary of the stable region in the parameter plane of the Mathieu equation. Furthermore, Kumar (1996) reduced the complexity of the mode significantly by adding the energy decay rate of the two-layer fluid to the Mathieu equations as a correction term. Labrador et al. (2021) derived the critical onset of Faraday waves, which appear via Hopf bifurcation, by extending the Mathieu equations to a three-layer system. However, this model was applied only to a specific case of two fluid layers alternating in microgravity, where all interfaces have the same natural frequencies and damping. The extension of this model to general cases may be of great significance for the understanding of Faraday waves considering viscous effects, which is one of the main objectives of the present study.

Since the discovery of Faraday waves, it has been accepted generally that the primary frequency of wave oscillation is half of the external excitation frequency, which is the basis of the main characteristic asymptotic analysis and Floquet analysis. Taking into account the effect of viscosity, Rajchenbach & Clamond (2015) added a viscosity correction term to the Mathieu equation and performed an asymptotic expansion and stability analysis. Their results suggest that the asymptotic analysis can provide a specific mathematical expression when estimating viscous dissipation as a small disturbance. However, it is important to note that this expression may contain significant errors when applied to high-viscosity fluids. The Mathieu equation can be solved using either asymptotic expansion or Floquet analysis. As a kind of numerical method, the Floquet analysis can reach a very high accuracy by using high-order Floquet expansion. It can solve the stable region of the viscous fluid in wavenumber space numerically. Kumar (1996) proposed this method to solve the complete fluid dynamics equations for hydrodynamical problems with single-layer or two-layer fluid. Concerning the liquid–liquid–gas system, Pototsky & Bestehorn (2016), Pototsky et al. (2004, 2019) and Bestehorn & Pototsky (2016) adopted the Floquet analysis method to study theoretically the dispersion relation and Rayleigh–Taylor instabilities. Recently, Ward, Zoueshtiagh & Narayanan (2019b) developed a complete three-layer viscous fluid model based on the Floquet analysis, and found double-tongued stability curves. However, the complexity of the model multiplies as the number of fluid layers increases. In addition, numerical methods can often encounter numerical singularities for fluid systems with extremely low viscosity under high-frequency vibrations when considering the complete fluid dynamics equations. Therefore, performing Floquet analysis on the Mathieu equation is a good choice to solve this problem. However, due to the unique nature of the three-layer fluid system with two fluid interfaces, it becomes necessary to solve directly the Mathieu equations numerically to ensure the reliability of analyses.

In the present study, we adopt asymptotic analysis, Floquet analysis and direct numerical simulation to study the Faraday instabilities of a three-layer fluid system. The purpose of the present paper has been threefold. First, we conduct the asymptotic expansion of Mathieu equations and Floquet analysis to obtain the properties of Faraday waves in a three-layer fluid system. Second, by solving Mathieu equations directly, we confirm the above analyses and identify the beating Faraday instability mode. Third, we clarify the occurrence of different modes and the connections between different modes.

We begin this paper by introducing the Mathieu equations governing Faraday waves of the three-layer system from inviscid to viscous fluids in § 2. The Floquet analysis of Mathieu equations is conducted in § 3, while the asymptotic expansion of Mathieu equations, the dispersion relation and the Faraday threshold are presented in § 4. In § 5, we build a Faraday experimental system and show the comparisons between experimental data and theoretical predictions. By solving the Mathieu equations directly, several Faraday instability modes of a double-interface system are discussed in § 6. Finally, we present some concluding remarks in § 7.

2. The Mathieu equations

We consider a cylindrical container filled with three Newtonian fluids of different densities, as depicted in figure 1(a). The motion of the three-layer fluid is governed by the Navier–Stokes equation

(2.1)\begin{equation} \rho_j \left (\frac{\partial \boldsymbol{V}_j}{t} + \boldsymbol{V}_j \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{V}_j \right)=-\boldsymbol{\nabla} P_j + \mu_j\,\nabla^2 \boldsymbol{V}_j-\rho_j \left[ g + a_c \cos 2\omega t \right] \boldsymbol{e}_z, \end{equation}

and the continuity equation

(2.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{V}_j = 0, \end{equation}

where $\boldsymbol {V}_j=u_j \boldsymbol {e}_x+v_j \boldsymbol {e}_y+w_j \boldsymbol {e}_z$, $P_j$, $\mu _j$ and $\rho _j$ are the velocity, pressure, kinetic viscosity and density of corresponding fluid $j$. Here, the subscript $j = 1, 2, 3$ denotes the liquid at the bottom, middle and top, respectively; $g$ denotes the acceleration due to gravity, while $a_c$ stands for the acceleration amplitude provided by the vibrator. In general, the vibration acceleration satisfies the formula $a_c = 4A \omega ^2$, where $A$ is the vibrational amplitude.

According to linear stability analysis, the solution of the ideal fluid equations (2.1)–(2.2) can be expanded in an asymptotic series as $\boldsymbol {V}_j=\boldsymbol {V}^{\prime }_j$, $P_j=P_{0,j}+P^{\prime }_j$. Therefore, the base state domain equations and the perturbation equations can be written as

(2.3)$$\begin{gather} \frac{\partial{P_{0,j}}}{\partial{z}}=-\rho_j \left(g+4A\omega^2 \cos\varOmega t\right), \end{gather}$$

(2.4)$$\begin{gather}\rho_j\,\frac{\partial{\boldsymbol{V}^{\prime}_{j}}}{\partial{t}}=-\boldsymbol{\nabla} P^{\prime}_j+\mu_j\,\nabla^2\boldsymbol{V}^{\prime}_j. \end{gather}$$

2.1. The inviscid model

Considering the ideal fluids, the non-rotating flow, and the stress-free sidewall assumption, similar to the Benjamin derivation process (Benjamin & Ursell 1954), we can get the $z$-direction velocity $w_j$ to satisfy the Helmholtz equation

(2.5)\begin{equation} \left( \frac{\partial^2}{\partial z^2}-k^2 \right) w_j=0, \end{equation}

where $k$ is the wavenumber in the horizontal plane. For cylindrical containers, $k=k_{ln}$ is the $n$th zero of ${\rm J}'_l(k_{ln}R)$, where ${\rm J}_l$ is a Bessel function of $l$th order. Equation (2.5) has a general solution

(2.6)\begin{equation} w_j=A_j(t)\mathrm{cosh}\,kz+B_j(t)\mathrm{sinh}\,kz. \end{equation}

At the wall, the no-penetration conditions are

(2.7)\begin{equation} \left.\begin{gathered} w_1=0 \quad \mathrm{at} \ z=0,\\ w_3=0 \quad \mathrm{at} \ z=h_1+h_2+h_3. \end{gathered}\right\}\end{equation}

At interface $A$, the flow satisfies the conditions

(2.8)$$\begin{gather} w_1=w_2=\frac{\mathrm{d}\xi_A}{\mathrm{d}t}, \end{gather}$$

(2.9)$$\begin{gather}\rho_2\,\frac{\partial}{\partial z}\frac{\partial w_2}{\partial t}-\rho_1\,\frac{\partial}{\partial z}\frac{\partial w_1}{\partial t} =-\left[ \left(\rho_2-\rho_1 \right)\left( g-a_c \cos 2\omega t\right)-\sigma_A k^2 \right]k^2 \xi_A. \end{gather}$$

At interface $B$, the flow satisfies the conditions

(2.10)$$\begin{gather} w_2=w_3=\frac{\mathrm{d}\xi_B}{\mathrm{d}t}, \end{gather}$$

(2.11)$$\begin{gather}\rho_3\,\frac{\partial}{\partial z}\frac{\partial w_3}{\partial t}-\rho_2\,\frac{\partial}{\partial z}\frac{\partial w_2}{\partial t} =-\left[ \left(\rho_3-\rho_2 \right)\left( g-a_c \cos 2\omega t\right)-\sigma_B k^2 \right]k^2 \xi_B. \end{gather}$$

From (2.6)–(2.11), the Mathieu equations of a three-layer inviscid fluid are obtained as

(2.12)$$\begin{gather} \frac{\mathrm{d}^2 \xi_A}{\mathrm{d} t^2}-\beta_A\,\frac{\mathrm{d}^2 \xi_B}{\mathrm{d} t^2} + \varOmega^2_A \left(1+F_A \cos 2\omega t\right) \xi_A = 0, \end{gather}$$

(2.13)$$\begin{gather}\frac{\mathrm{d}^2 \xi_B}{\mathrm{d} t^2}-\beta_B\,\frac{\mathrm{d}^2 \xi_A}{\mathrm{d} t^2} + \varOmega^2_B \left(1+F_B \cos 2\omega t\right) \xi_B = 0, \end{gather}$$

where $\varOmega _A$ and $\varOmega _B$ are natural frequencies of the interfaces $A$ and $B$:

(2.14a,b)\begin{equation} \varOmega^2_A=\frac{(\rho_1-\rho_2)gk+\sigma_A k^3}{\rho_1\,\mathrm{coth}\, kh_1+\rho_2\,\mathrm{coth}\,kh_2},\quad \varOmega^2_B=\frac{(\rho_2-\rho_3)gk+\sigma_B k^3}{\rho_2\,\mathrm{coth} \,kh_2+\rho_3\,\mathrm{coth}\,kh_3}. \end{equation}

In addition, other parameters are

(2.15a,b)$$\begin{gather} F_A = \frac{\left(\rho_1-\rho_2\right)a_c}{\sigma_A k^2+\left(\rho_1-\rho_2\right)g},\quad F_B = \frac{\left(\rho_2-\rho_3\right)a_c}{\sigma_B k^2+\left(\rho_2-\rho_3\right)g}, \end{gather}$$

(2.16a,b)$$\begin{gather}\beta_A = \frac{\rho_2/\mathrm{sinh}\,kh_2}{\rho_1\,\mathrm{coth}\,kh_1+\rho_2\,\mathrm{coth}\,kh_2},\quad \beta_B = \frac{\rho_2/\mathrm{sinh}\,kh_2}{\rho_2\,\mathrm{coth}\,kh_2+\rho_3\,\mathrm{coth}\,kh_3}, \end{gather}$$

where subscripts $A$ and $B$ represent the interfaces $A$ and $B$, respectively. Here, $F_A$ and $F_B$ describe the relative magnitudes of vibration acceleration, and $\beta _A$ and $\beta _B$ describe the strength of coupling between interfaces.

We then derive the Mathieu equations of an $N$-layer fluid model (figure 2) for the inviscid fluids by adopting a method similar to that for the three-layer fluid model, and present the details in Appendix A.

Figure 2. Schematic diagram of $N$-layer fluid system.

The following simplified equations (2.17), (2.18) can be obtained easily from (2.12), (2.13) under the assumption that $kh_2\gg 1$:

(2.17)$$\begin{gather} \frac{\mathrm{d}^2\xi_A}{\mathrm{d}t^2}+\varOmega^2_A\left( 1+F_A \cos 2\omega t \right)\xi_A=0, \end{gather}$$

(2.18)$$\begin{gather}\frac{\mathrm{d}^2\xi_B}{\mathrm{d}t^2}+\varOmega^2_B\left( 1+F_B \cos 2\omega t \right)\xi_B=0. \end{gather}$$

In this case, the equations are decoupled into two two-layer inviscid models, indicating that there is no coupling between the motions of the two interfaces when the thickness of the medium-density fluid is large enough. It is consistent with the results from Ward et al. (2019b).

2.2. The viscous model

While inviscid models can explain partially the mechanism of Faraday waves, they cannot provide the acceleration threshold due to the neglect of viscous dissipation. Landau & Lifshitz (1987) presented a method to calculate the rate of viscous dissipation by estimating the energy of viscous dissipation. The Faraday wave of a two-layer fluid is extended further and compared with Floquet's model (Kumar & Tuckerman 1994), showing good agreement. This study proposes a simplification in a three-layer viscous model (Labrador et al. 2021). The Mathieu equations of a three-layer viscous fluid can be obtained by considering only the direct viscous action of fluids on both sides of the interface:

(2.19)$$\begin{gather} \frac{\mathrm{d}^2\xi_A}{\mathrm{d}t^2}-\beta_A\,\frac{\mathrm{d}^2\xi_B}{\mathrm{d}t^2}+2\gamma_A\, \frac{\mathrm{d}\xi_A}{\mathrm{d}t}+\varOmega^2_A\left( 1+F_A\cos 2\omega t\right)\xi_A=0, \end{gather}$$

(2.20)$$\begin{gather}\frac{\mathrm{d}^2\xi_B}{\mathrm{d}t^2}-\beta_B\,\frac{\mathrm{d}^2\xi_A}{\mathrm{d}t^2}+2\gamma_B\, \frac{\mathrm{d}\xi_B}{\mathrm{d}t}+\varOmega^2_B\left( 1+F_B\cos 2\omega t\right)\xi_B=0, \end{gather}$$

where $\gamma _A$ and $\gamma _B$ are the damping coefficients of the interfaces $A$ and $B$,

(2.21a,b)\begin{align} \gamma_A = \frac{2k^2\left( \mu_1\,\mathrm{coth}\,kh_1+\mu_2 \,\mathrm{coth}\,kh_2\right)}{\rho_1\,\mathrm{coth}\,kh_1+\rho_2\,\mathrm{coth}\,kh_2},\quad \gamma_B = \frac{2k^2\left( \mu_2\,\mathrm{coth}\,kh_2+\mu_3 \,\mathrm{coth}\,kh_3\right)}{\rho_2\,\mathrm{coth}\,kh_2+\rho_3\,\mathrm{coth}\,kh_3}. \end{align}

Applying the Faraday instability in microgravity, the viscous model ((2.19) and (2.20)) is validated by the simulations from Labrador et al. (2021). We can estimate these damping coefficients through an evaluation of viscous dissipation energy. The longer the wavelength, the more significant the viscous damping (Batson et al. 2013) of the interface and the wall boundary layer becomes. When the wavenumber is large, viscous damping arises primarily from internal damping, which can be estimated by calculating the dissipated viscous energy (Herreman et al. 2019; Kumar & Tuckerman 1994).

On the basis of (A7), we add the viscous correction term to get the Mathieu equations of an $n$-layer viscous fluid:

(2.22)\begin{equation} \begin{cases} \dfrac{\text{d}^2\xi_1}{\text{d}t^2}-\beta_{1,2}\,\dfrac{\text{d}^2\xi_2}{\text{d}t^2}+2\gamma_1\, \dfrac{\mathrm{d}\xi_1}{\mathrm{d}t}+\varOmega_1\left( 1+F_1\cos 2\omega t \right)\xi_1 =0,\\[9pt] \vdots\\[4pt] \dfrac{\text{d}^2\xi_i}{\text{d}t^2}-\beta_{i,i-1}\,\dfrac{\text{d}^2\xi_{i-1}}{\text{d}t^2}-\rho_{i,i+1}\,\dfrac{\text{d}^2\xi_{i+1}}{\text{d}t^2}+2\gamma_i\, \dfrac{\mathrm{d}\xi_i}{\mathrm{d}t}+\varOmega_i\left( 1+F_i\cos 2\omega t \right)\xi_i =0, \\[9pt] \quad 1< i< N-1,\\[3pt] \vdots\\[4pt] \dfrac{\text{d}^2\xi_{N-1}}{\text{d}t^2}-\beta_{N-1,N-2}\,\dfrac{\text{d}^2\xi_{N-2}}{\text{d}t^2}+2\gamma_{N-1}\, \dfrac{\mathrm{d}\xi_{N-1}}{\mathrm{d}t}+\varOmega_{N-1}\left( 1+F_{N-1}\cos 2\omega t \right)\xi_{N-1} =0, \end{cases} \end{equation}

where

(2.23)\begin{equation} \gamma_i=\frac{2k^2\left( \mu_i \,\mathrm{coth}\,kh_i+\mu_{i+1} \,\mathrm{coth}\,kh_{i+1} \right)}{\rho_i \,\mathrm{coth}\,kh_i+\rho_{i+1} \,\mathrm{coth}\,kh_{i+1}}. \end{equation}

3. Floquet analysis

Floquet analysis can be applied to parametric instability problems such as Faraday waves. Ward et al. (2019b) performed a Floquet analysis of the three-layer Faraday waves on the Kumar theoretical model (Kumar & Tuckerman 1994; Kumar 1996), and found that Floquet analysis can be generalized to three-layer fluid systems. In the same way, we carried out Floquet analysis on the Mathieu equations (2.19) of a three-layer viscous fluid.

The Floquet expansion is used to separate the variables

(3.1)\begin{equation} \left.\begin{gathered} \xi_A = {\rm e}^{{\rm i} kx} \sum_{n=-\infty }^{\infty } {\rm e}^{q_n t}\,\hat{\xi} _{A,n},\\ \xi_B = {\rm e}^{{\rm i} kx} \sum_{n=-\infty }^{\infty } {\rm e}^{q_n t}\,\hat{\xi} _{B,n}, \end{gathered}\right\} \end{equation}

where $q_n=s+\mathrm {i}( \alpha +n)2\omega$. The stability curves for harmonic ($\alpha = 0$) and subharmonic ($\alpha = 1/2$) solutions can be obtained by setting $s=0$, where $s$ is the growth rate. Substitute (3.1) into (2.19)–(2.20) to get

(3.2)\begin{equation} \left.\begin{gathered} \left( q^2_n+2q_n \gamma_A+\varOmega^2_A\right)\hat{\xi}_{A,n}-\beta_A q^2_n\, \hat{\xi}_{B,n}+\frac{\varOmega^2_A F_A}{2} \left( \hat{\xi}_{A,n-1}+\hat{\xi}_{A,n+1}\right)=0,\\ \left( q^2_n+2q_n \gamma_B+\varOmega^2_B\right)\hat{\xi}_{B,n}-\beta_B q^2_n\, \hat{\xi}_{A,n}+\frac{\varOmega^2_B F_B}{2} \left( \hat{\xi}_{B,n-1}+\hat{\xi}_{B,n+1}\right)=0. \end{gathered}\right\} \end{equation}

These equations can be formulated as an eigenvalue problem for $\hat {\xi }_{A,n}, \hat {\xi }_{B,n}$, where $A$ is the eigenvalue. The eigenvalue equation is

(3.3)\begin{equation} \begin{pmatrix} q^2_n+2q_n\gamma_A+\varOmega^2_A & -\beta_A q^2_n \\ -\beta_B q^2_n & q^2_n+2q_n\gamma_B+\varOmega^2_B \end{pmatrix} \begin{pmatrix} \hat{\xi}_{A,n} \\ \hat{\xi}_{B,n} \end{pmatrix} = A\omega^2 \begin{pmatrix} \chi_A & 0 \\ 0 & \chi_B \end{pmatrix} \begin{pmatrix} \hat{\xi}_{A,n-1}+\hat{\xi}_{A,n+1} \\ \hat{\xi}_{B,n-1}+\hat{\xi}_{B,n+1} \end{pmatrix}\!, \end{equation}

where

(3.4a,b)\begin{equation} \chi_A =-\frac{\varOmega^2_A}{2}\,\frac{\rho_1-\rho_2}{\sigma_A k^2+( \rho_1-\rho_2)g},\quad \chi_B =-\frac{\varOmega^2_B}{2}\,\frac{\rho_2-\rho_3}{\sigma_B k^2+( \rho_2-\rho_3)g}. \end{equation}

Further, we can write (3.3) as a system of linear equations

(3.5)\begin{equation} \boldsymbol{M} \begin{pmatrix} \hat{\xi}_A \\ \hat{\xi}_B \end{pmatrix} = A \omega^2 \begin{pmatrix} \chi_A \boldsymbol{C} & 0 \\ 0 & \chi_B \boldsymbol{C} \end{pmatrix} \begin{pmatrix} \hat{\xi}_A \\ \hat{\xi}_B \end{pmatrix}, \end{equation}

where $\hat {\xi }_A,\hat {\xi }_B$ are column vectors

(3.6a,b)\begin{equation} \hat{\xi}_A=\begin{pmatrix} \hat{\xi}_{A,-n} \\ \hat{\xi}_{A,-n+1} \\ \vdots \\ \hat{\xi}_{A,0} \\ \vdots \\ \hat{\xi}_{A,n-1} \\ \hat{\xi}_{A,n} \end{pmatrix},\quad \hat{\xi}_B=\begin{pmatrix} \hat{\xi}_{B,-n} \\ \hat{\xi}_{B,-n+1} \\ \vdots \\ \hat{\xi}_{B,0} \\ \vdots \\ \hat{\xi}_{B,n-1} \\ \hat{\xi}_{B,n} \end{pmatrix}, \end{equation}

and $\boldsymbol {M}$ and $\boldsymbol {C}$ are banded matrices

(3.7)\begin{gather} \boldsymbol{M} = \begin{pmatrix}\begin{smallmatrix} q^2_{-n}+2q_{-n}\gamma_A+\varOmega^2_A & & & & & & -\beta_A q^2_{-n} & & & & & \\ & & & \ddots & & & & & & \ddots & & \\ & & & & & q^2_n+2q_n\gamma_A+\varOmega^2_A & & & & & & -\beta_A q^2_n\\ -\beta_B q^2_{-n} & & & & & & q^2_{-n}+2q_{-n}\gamma_B+\varOmega^2_B & & & & & \\ & & & \ddots & & & & & & \ddots & & \\ & & & & & -\beta_B q^2_n & & & & & & q^2_n+2q_n\gamma_B+\varOmega^2_B \end{smallmatrix}\end{pmatrix}, \end{gather}

(3.8)\begin{gather} \boldsymbol{C} = \begin{pmatrix} 0 & 1 & & & & \\ 1 & 0 & 1 & & & \\ & 1 & 0 & 1 & & \\ & & 1 & 0 & 1 & \\ & & & 1 & 0 & 1 \\ & & & & 1 & 0 \end{pmatrix}. \end{gather}

The stability boundaries are defined by the curves in the $A\unicode{x2013} k$ plane on which $s(A, k)=0$, for a certain frequency $\omega$. Note that the eigenvalue $A$ must be a positive real number.

4. Asymptotic analysis

To facilitate the solution, the dimensionless Mathieu equations of a three-layer viscous fluid (2.19)–(2.20) are obtained, where the time scale is $1/\omega$:

(4.1)$$\begin{gather} \ddot{\xi}_A - \beta_A \ddot{\xi}_B + 2\hat{\mu}_A \dot{\xi}_A + \left( \delta_A + \hat{\varepsilon}_A \cos 2t \right)\xi_A=0, \end{gather}$$

(4.2)$$\begin{gather}\ddot{\xi}_B - \beta_B \ddot{\xi}_A + 2\hat{\mu}_B \dot{\xi}_B + \left( \delta_B + \hat{\varepsilon}_B \cos 2t \right)\xi_B=0, \end{gather}$$

where $\hat {\mu }_A = \gamma _A/\omega$, $\hat {\mu }_B = \gamma _B/\omega$, $\delta _A = \varOmega ^2_A/\omega ^2$, $\delta _B = \varOmega ^2_B/\omega ^2$, $\hat {\varepsilon }_A = F_A\varOmega ^2_A/\omega ^2$ and $\hat {\varepsilon }_B = F_B\varOmega ^2_B/\omega ^2$.

Let a small quantity of the system be $\epsilon = a_c/g$. It is also assumed that $\hat {\mu }_A$ and $\hat {\mu }_B$ are of order $O ( \epsilon )$ for a low-viscosity fluid system. Then $\hat {\mu }_A=\epsilon \mu _A$, $\hat {\mu }_B=\epsilon \mu _B$, $\hat {\varepsilon }_A=\epsilon \varepsilon _A$ and $\hat {\varepsilon }_B=\epsilon \varepsilon _B$. Physically, $\epsilon$ is the dimensionless amplitude of externally excited acceleration. This treatment method is based on the work of Rajchenbach & Clamond (2015). Asymptotically expand $\xi _A$, $\xi _B$, $\delta _A$ and $\delta _B$, as

(4.3)\begin{equation} \left.\begin{gathered} \xi_A=\xi_{A0}+\epsilon \xi_{A1}+\epsilon^2 \xi_{A2}+O (\epsilon^3), \\ \xi_B=\xi_{B0}+\epsilon \xi_{B1}+\epsilon^2 \xi_{B2}+O (\epsilon^3), \\ \delta_A=\delta_{A0}+\epsilon \delta_{A1}+\epsilon^2 \delta_{A2}+O (\epsilon^3), \\ \delta_B=\delta_{B0}+\epsilon \delta_{B1}+\epsilon^2 \delta_{B2}+O (\epsilon^3). \end{gathered}\right\} \end{equation}

According to the characteristics of Faraday waves, at order $O (\epsilon ^0)$, the natural frequency of the fluid system is about half of the excitation frequency, that is, $\omega \approx \omega _0$, where $\omega _0$ is the natural frequency of the fluid system. Therefore, we define $\delta _{A0} = \varOmega _A ^ 2 / \omega _0 ^ 2$ and $\delta _{B0} = \varOmega _B ^ 2 / \omega _0 ^ 2$.

By substituting (4.3) into (4.1)–(4.2), it is easy to get the zero-order approximation equations (4.4), the first-order approximation equations (4.5), and the second-order approximation equations (4.6):

(4.4)\begin{gather} \left.\begin{gathered} \ddot{\xi}_{A0}-\beta_A \ddot{\xi}_{B0} +\delta_{A0} \xi_{A0} = 0,\\ \ddot{\xi}_{B0}-\beta_B \ddot{\xi}_{A0} +\delta_{B0} \xi_{B0} = 0.\end{gathered}\right\} \end{gather}

(4.5)\begin{gather} \left.\begin{gathered} \ddot{\xi}_{A1}-\beta_A \ddot{\xi}_{B1} +\delta_{A0} \xi_{A1} =-2\mu_A \dot{\xi}_{A0}-\left( \delta_{A1}+ \varepsilon_A \cos 2t \right)\xi_{A0},\\ \ddot{\xi}_{B1}-\beta_B \ddot{\xi}_{A1} +\delta_{B0} \xi_{B1} =-2\mu_B \dot{\xi}_{B0}-\left( \delta_{B1}+ \varepsilon_B \cos 2t \right)\xi_{B0}. \end{gathered}\right\} \end{gather}

(4.6)\begin{gather} \left.\begin{gathered} \ddot{\xi}_{A2}-\beta_A \ddot{\xi}_{B2} +\delta_{A0} \xi_{A2} =-2\mu_A \dot{\xi}_{A1}-\delta_{A2}\xi_{A0}- \left( \delta_{A1}+ \varepsilon_A \cos 2t \right)\xi_{A1},\\ \ddot{\xi}_{B2}-\beta_B \ddot{\xi}_{A2} +\delta_{B0} \xi_{B2} =-2\mu_B \dot{\xi}_{B1}-\delta_{B2}\xi_{B0}-\left( \delta_{B1}+ \varepsilon_B \cos 2t \right)\xi_{B1}. \end{gathered}\right\} \end{gather}

The solution of (4.4) can be written in the form

(4.7)\begin{equation} \left.\begin{gathered} \xi_{A0} = A_0\cos t+B_0\sin t,\\ \xi_{B0} = C_0\cos t+D_0\sin t, \end{gathered}\right\} \end{equation}

where $A_0$, $B_0$, $C_0$, $D_0$ are all undetermined coefficients. Substitute (4.7) into (4.4) to derive the subsequent equations

(4.8a)$$\begin{gather} \begin{pmatrix} \delta_{A0}-1 & \beta_A\\ \beta_B & \delta_{B0}-1 \end{pmatrix} \begin{pmatrix}{A_0}\\ {C_0}\end{pmatrix} =0, \end{gather}$$

(4.8b)$$\begin{gather}\begin{pmatrix} \delta_{A0}-1 & \beta_A\\ \beta_B & \delta_{B0}-1 \end{pmatrix} \begin{pmatrix}{B_0}\\ {D_0}\end{pmatrix} =0. \end{gather}$$

Since the undetermined coefficients ($A_0, B_0, C_0, D_0$) cannot be equal to 0, which implies the existence of a non-zero solution, the determinant of the coefficient matrix must be 0:

(4.9)\begin{equation} \begin{vmatrix} \delta_{A0}-1 & \beta_A \\ \beta_B & \delta_{B0}-1 \end{vmatrix} =0. \end{equation}

Equation (4.9) can also be written as

(4.10)\begin{equation} \begin{vmatrix} \dfrac{\varOmega_A ^ 2}{\omega_0^ 2} -1 & \beta_A \\[10pt] \beta_B & \dfrac{\varOmega_B ^ 2}{\omega_0^2}-1 \end{vmatrix} =0. \end{equation}

In this case, $\omega _0$ satisfies

(4.11)\begin{equation} \beta_A\beta_B=\frac{\left( \omega^2_0-\varOmega^2_A \right)\left( \omega^2_0-\varOmega^2_B \right)}{\omega^4_0}, \end{equation}

where $\omega _{0,1}$ and $\omega _{0,2}$ ($\omega _{0,1}>\omega _{0,2}$) are the two non-negative roots of the equation.

Substitute (4.7) into (4.5) to get

(4.12) \begin{align} \ddot{\xi}_{A1}-\beta_A \ddot{\xi}_{B1}+\delta_{A0}\xi_{A1} &= \frac{1}{2}\left(-2A_0\delta_{A1}-A_0\varepsilon_A-4B_0\mu_A\right)\cos t\nonumber\\ &\quad +\frac{1}{2}\left(-2B_0\delta_{A1}+B_0\varepsilon_A+4A_0\mu_A \right)\sin t\nonumber\\ &\quad +\frac{1}{2}\left(-A_0 \varepsilon_A \cos 3t-B_0\varepsilon_A\sin 3t \right), \end{align}

(4.13)\begin{align} \ddot{\xi}_{B1}-\beta_B \ddot{\xi}_{A1}+\delta_{B0}\xi_{B1} &= \frac{1}{2}\left(-2C_0\delta_{B1}-C_0\varepsilon_B-4D_0\mu_B\right)\cos t \nonumber\\ &\quad +\frac{1}{2}\left(-2D_0\delta_{B1}+D_0\varepsilon_B+4C_0\mu_B \right)\sin t\nonumber\\ &\quad +\frac{1}{2}\left(-C_0 \varepsilon_B \cos 3t-D_0\varepsilon_B\sin 3t \right). \end{align}

The secular terms in (4.12) and (4.13) are eliminated to obtain

(4.14a,b)\begin{equation} \delta_{A1}={\pm}\sqrt{\frac{\varepsilon^2_A}{4}-4\mu^2_A},\quad \delta_{B1}={\pm}\sqrt{\frac{\varepsilon^2_B}{4}-4\mu^2_B}, \end{equation}

and

(4.15)\begin{equation} \left.\begin{gathered} \xi_{A1}=\frac{A_0\varepsilon_A\left(\delta_{B0}-9\right)-9C_0\beta_A\varepsilon_B}{162 \beta_A \beta_B-2\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)}\cos 3t + \frac{B_0\varepsilon_A\left(\delta_{B0}-9\right)-9D_0\beta_A\varepsilon_B}{162 \beta_A \beta_B-2\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)}\sin 3t,\\ \xi_{B1}=\frac{C_0\varepsilon_B\left(\delta_{A0}-9\right)-9A_0\beta_B\varepsilon_A}{162 \beta_A \beta_B-2\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)}\cos 3t + \frac{D_0\varepsilon_B\left(\delta_{A0}-9\right)-9B_0\beta_B\varepsilon_B}{162 \beta_A \beta_B-2\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)}\sin 3t. \end{gathered}\right\} \end{equation}

In the same way, the secular terms of (4.6) are eliminated to obtain

(4.16)\begin{equation} \left.\begin{gathered} \delta_{A2} = \frac{\varepsilon_A\left(\left(\delta_{A0}-1\right)\left(\delta_{B0}-9\right)\varepsilon_A+9\beta_A^2\varepsilon_B\right)}{4\left(\delta_{A0}-1\right)\left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)-81\beta_A\beta_B\right)},\\ \delta_{B2} = \frac{\varepsilon_B\left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-1\right)\varepsilon_B+9\beta_B^2\varepsilon_A\right)}{4\left(\delta_{B0}-1\right)\left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)-81\beta_A\beta_B\right)}. \end{gathered}\right\} \end{equation}

From (4.14a,b) and (4.16), the second-order approximate solutions for $\delta _A$ and $\delta _B$ are obtained as

(4.17)\begin{equation} \left.\begin{gathered} \delta_A = \frac{\varOmega^2_A}{\omega^2_0}\pm\sqrt{\frac{\hat{\varepsilon}^2_A}{4}-4\mu^2_A}+ \frac{\hat{\varepsilon}_A\left(\left(\delta_{A0}-1\right)\left(\delta_{B0}-9\right) \hat{\varepsilon}_A+9\beta_A^2\hat{\varepsilon}_B\right)}{4\left(\delta_{A0}-1\right) \left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)-81\beta_A\beta_B\right)},\\ \delta_B = \frac{\varOmega^2_B}{\omega^2_0}\pm\sqrt{\frac{\hat{\varepsilon}^2_B}{4}- 4\mu^2_B}+\frac{\hat{\varepsilon}_B\left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-1\right) \hat{\varepsilon}_B+9\beta_B^2\hat{\varepsilon}_A\right)}{4\left(\delta_{B0}-1\right) \left(\left(\delta_{A0}-9\right)\left(\delta_{B0}-9\right)-81\beta_A\beta_B\right)}. \end{gathered}\right\} \end{equation}

It can be seen that the interface coupling effect exists not only in the dispersion relation (4.11), but also in the second-order term of acceleration (4.17). However, viscosity is one of the factors to inhibit Faraday instability at both multi-fluid and single-fluid interfaces.

Furthermore, (4.14a,b) can be written in terms of variations of $\omega$:

(4.18)\begin{equation} \left( 1-\frac{\omega^2}{\omega_0^2} \right) \approx \frac{F_A^2}{4} - \frac{4\gamma_A^2 \omega^2}{\varOmega_A^4}=O\left( \epsilon \right) \end{equation}

and

(4.19)\begin{equation} \left( 1-\frac{\omega^2}{\omega_0^2} \right) \approx \frac{F_B^2}{4} - \frac{4\gamma_B^2 \omega^2}{\varOmega_B^4}=O\left( \epsilon \right). \end{equation}

So when $\omega \approx \omega _0$ is satisfied, critical instability occurs. From (4.14a,b), we can obtain

(4.20a,b)\begin{equation} F_A=\frac{4\gamma_A}{\varOmega^2_A}\omega_0,\quad F_B=\frac{4\gamma_B}{\varOmega^2_B}\omega_0. \end{equation}

Further, (4.20a,b) is dimensionalized to get

(4.21)\begin{equation} \left.\begin{gathered} a_{cA}=8k\omega\,\frac{\mu_1\,\mathrm{coth}\,kh_1+\mu_2\,\mathrm{coth}\,kh_2}{\rho_1-\rho_2},\\ a_{cB}=8k\omega\,\frac{\mu_2\,\mathrm{coth}\,kh_2+\mu_3\,\mathrm{coth}\,kh_3}{\rho_2-\rho_3}, \end{gathered}\right\} \end{equation}

where $a_{cA}$ and $a_{cB}$ correspond to the critical excitation accelerations of two subharmonic instability modes (discussed in detail in § 6), respectively. To facilitate the comparison with the experimental results, the vibration amplitude (4.22) is derived:

(4.22a)$$\begin{gather} A_{cA}=\frac{2k}{\omega}\,\frac{\mu_1\,\mathrm{coth}\,kh_1+\mu_2\,\mathrm{coth}\,kh_2}{\rho_1-\rho_2}, \end{gather}$$

(4.22b)$$\begin{gather}A_{cB}=\frac{2k}{\omega}\,\frac{\mu_2\,\mathrm{coth}\,kh_2+\mu_3\,\mathrm{coth}\,kh_3}{\rho_2-\rho_3} . \end{gather}$$

It is worth noting that when $kh_2\gg 1$ (i.e. the middle fluid is thick enough), a result identical to that of the single interface is obtained:

(4.23a,b)\begin{equation} \omega_{0,A}=\varOmega_A,\quad \omega_{0,B}=\varOmega_B \end{equation}

and

(4.24)\begin{equation} \left.\begin{gathered} \delta_A = \frac{\varOmega^2_A}{\omega^2_0}\pm\sqrt{\frac{\hat{\varepsilon}^2_A}{4}-4\mu^2_A}+\frac{\hat{\varepsilon}_A^2}{4 (\delta_{A0}-9)}, \\ \delta_B = \frac{\varOmega^2_B}{\omega^2_0}\pm\sqrt{\frac{\hat{\varepsilon}^2_B}{4}-4\mu^2_B}+\frac{\hat{\varepsilon}_B^2}{4 (\delta_{B0}-9)}. \end{gathered}\right\} \end{equation}

Equations (4.23a,b)–(4.24) are the solution for interfaces $A$ and $B$. In this case, the double interface is completely decoupled, and the fluctuations do not interfere with each other.

5. Experimental systems and measurement methods

The schematic diagram of the double-interface fluid system is shown in figure 3. A cylindrical tank of inner diameter 55 mm is mounted on a vertical vibration system that consists of an SA-JZ020 electromagnetic shaker, an FY6900 function generator and an HEA-500G power amplifier. With the help of the CDAQ-9171 data acquisition system, the acceleration of the tank measured by the SAE D0005B accelerometer can be stored with a high sample rate. A homemade Python code is adopted to realize a combined operation of the vertical vibration system, the data acquisition system and the high-speed camera system (Phantom VEO 340S).

Figure 3. Schematic diagram of the experimental system.

We select eutectic alloy GaInSn (66.7 % Ga, 20.5 % In, 12.5 % Sn), which is in a liquid state at room temperature, as the working fluid. We place 30 ml (${h_1=12.62\ \textrm {mm}}$) GaInSn in the cylindrical container, and then add 3 ml (${h_2=1.262\ \textrm {mm}}$) sodium hydroxide solution ($1\ \textrm {mol}\ \textrm {l}^{-1}$) on the surface of the GaInSn to prevent oxidizing. Here, it is important to note that although the oxide surface layer can be removed by both acidic and alkaline solutions, we have chosen to use a $1\ \textrm {mol}\ \textrm {l}^{-1}$ concentration of sodium hydroxide due to its milder chemical reaction (Handschuh-Wang et al. 2018). The interfacial tension coefficients $\sigma _A$ between the GaInSn and the sodium hydroxide solution, and $\sigma _B$ between the sodium hydroxide solution and air, are obtained by the pendant drop method. The OCA25 contact angle measuring instrument (manufactured in Germany) measures the interfacial tension coefficient with resolution $0.01\ \textrm {mN}\ \textrm {m}^{-1}$. Since the air density is much smaller than in the GaInSn and sodium hydroxide solution, we consider only the physical properties of the GaInSn and sodium hydroxide solution that are listed in table 1.

Table 1. The fluid properties.

A high-speed camera is mounted directly above the container to capture the shape of the interface with the help of LED lights on the side of the container through a semi-transparent and semi-reflective mirror. Regarding the two fluid interfaces present in the system, it is difficult to distinguish each interface accurately when the fluids are transparent, as in Ward et al. (2019b) while studying the Faraday wave in three-layer fluid systems. In the present study, however, the situation is quite different. The reflectance of the GaInSn/solution interface is much stronger than that of the solution/air interface, and as a result, the camera captures most of the reflected light from the GaInSn/solution interface and less from the solution/air interface. Since the camera's optical axis is parallel to the reflected light, the bright areas correspond to the antinodes of the wave, while the dark areas correspond to the nodes. The magnitude of the brightness is determined by the slopes of the surface waves, with darker regions having larger slopes as shown in figure 4. This allows the modal pattern of the GaInSn/solution interface surface waves to be acquired by overhead recording, reflecting the spatial distribution of the surface waves. Moreover, in the present experimental system, two interfaces are coupled strongly to each other and exhibit synchronous fluctuations.

Figure 4. (a–c) Theoretical prediction of surface shapes of mode $(n, l)$. (d–f) Surface shapes plotted by one minus absolute value. (g–i) Faraday patterns. The driving frequency of the vibrator is $f$, which satisfies $\varOmega =2{\rm \pi} f$, where $n$ and $l$ are the radial and azimuthal mode numbers.

Due to the wetting properties, a meniscus is formed at the interface between the GaInSn and the sodium hydroxide solution, curving towards the GaInSn side. Similarly, at the interface between the sodium hydroxide solution and air, a meniscus is formed curving towards the air side. As suggested by Shao et al. (2021a,b,c) and Wilson et al. (2022), these menisci produce axisymmetric and harmonic edge waves in the experiment that change only the modes with $l=0$. In this study, we focus on situations where $l > 1$, therefore we can ignore the edge waves generated by the contact angle.

To achieve a high-resolution wave mode, we scanned with 0.1 Hz frequency interval and 0.1 $\mathrm {m}\ \mathrm {s}^{-2}$ acceleration interval during experiments. Within a certain range of frequency and acceleration, a type of Faraday wave occurs. In experiments, the characteristic excitation frequency is defined as the excitation frequency corresponding to the minimum acceleration, which is approximately where the edge wave dominates and the interface fluctuates slightly at low excitation accelerations below the critical value. However, when the excitation acceleration exceeds the critical value, the Faraday wave is generated and the amplitude of the interface wave increases significantly. It is worth noting that the Faraday wave transitions gradually to chaos as the excitation acceleration increases.

Since the boundary curve of the container is a circle, the spatial structure of Faraday waves (figures 4a–c) in a cylindrical container can be described as $\textrm {J}_l( k_{ln} r )\exp (\mathrm {i} l\theta )$ (Benjamin & Ursell 1954) due to the Helmholtz equation in polar form, i.e.

(5.1)\begin{equation} \left(\frac{\partial^2 }{\partial r^2}+\frac{1}{r}\,\frac{\partial }{\partial r}+\frac{1}{r^2}\, \frac{\partial^2 }{\partial \theta^2} + k_{ln}^2 \right)\xi_j=0 ,\quad j=A,B. \end{equation}

Here, for cylindrical containers, $k=k_{ln}$ is the $n$th zero of $\textrm {J}'_l(k_{ln}R)$, where $\textrm {J}_l$ is a Bessel function of $l$th order.

To determine the mode that corresponds to the observed patterns (figures 4g–i), we calculate one minus the absolute value of the wave slope, which is obtained by superimposing the images of the two periods of the Faraday waves, as shown in figures 4(d–f). In addition, the modal structure is identified by two-dimensional cross-correlation with the observed patterns (Shao et al. 2021c).

To determine the relative deformations of the two interfaces, the positions of the camera and lights are changed so that the contact lines of the two surfaces can be observed clearly, as shown in figure 5(e). We select four typical cases with different driving frequencies at which the Faraday wave mode can be observed. From figures 5(a–d), it can be seen that the waveforms of two interfaces comply with the definition of the in-phase deformation, corresponding to the definition of the zigzag mode in space. We further calculate the evolution of the heights of positions $A$ and $B$ marked in figure 5(d), and plot the time-dependent data in figure 5(f). It can be seen that the temporal evolutions of the two interfaces are synchronized. Therefore, we can conclude that the two interfaces fluctuate synchronously and exhibit the characteristics of the zigzag mode in experiments. Since two interfaces show the same pattern, we can obtain the interface pattern from either interface in a strongly coupled system.

Figure 5. (a–d) Side-view snapshots of four typical Faraday wave modes. The positions of interfaces $A$ and $B$ are marked by solid blue and red lines, respectively. (e) Schematic diagram of the arrangement of camera and container. (f) Temporal evolution of the Faraday wave amplitude at 62.6 Hz.

According to the asymptotic analysis, the dispersion relation of Faraday can be described approximately by (4.11), i.e. $\omega \approx \omega _0$. As shown in figure 6, for the zigzag mode, the experimental results and theoretical calculations from (4.11) agree well. Furthermore, the frequency increases as the wavenumbers ($l$ or $n$) increase, which is in line with previous research on Faraday waves (Benjamin & Ursell 1954; Rajchenbach & Clamond 2015; Pototsky & Bestehorn 2016; Shao et al. 2021c). However, there is a difference between experimental results and predicted data from the dispersion equation (4.11) regarding the varicose mode, which will be explained in detail in § 6.1 when the coupling effects of two interfaces are discussed.

Figure 6. Comparison of dispersion relations. The red dots are the experimental results of § 5. The blue lines are the calculated results from the dispersion equation of zigzag mode (4.11), while the green lines are from the dispersion equation of the varicose mode (4.11).

In addition, the vibration amplitude required for Faraday waves was acquired, ranging from 30 to 100 Hz with interval 0.2 Hz, as shown in figure 7. The given data points represent the minimum oscillation amplitude required to observe Faraday waves at a specific frequency. Due to the geometric constraints of the experimental vessel, the Faraday waves have discrete wavenumbers, which results in critical amplitudes of oscillation that are distributed over several Faraday tongues (critical stability curves). In figure 7, we plot the fitted three Faraday tongues corresponding to modes $(2,6)$, $(5,2)$, and $(6,7)$, respectively. The extremum points on each Faraday tongue have the minimum critical oscillation amplitude that corresponds to a particular mode. Figure 7 illustrates that as the frequency decreases, the Faraday tongue narrows, and the variation in critical oscillation amplitude becomes more pronounced. Additionally, the distance between Faraday tongues becomes smaller as the frequency increases, because the frequency differences between the modes become smaller. Therefore, a large amount of data is required to find exactly the minimum value of the Faraday tongue that corresponds to a mode.

Figure 7. Comparison of our experiments with the critical vibration amplitude-driven zigzag mode Faraday wave in Floquet analysis and asymptotic analysis. The dashed lines are the fitted Faraday tongues.

The corresponding critical vibration amplitudes driving the Faraday waves (Faraday threshold), from experiments, Floquet analysis on the Mathieu equations, first-order asymptotic analysis and second-order asymptotic analysis, are displayed in figure 7. One can see that the critical vibration amplitudes obtained by Floquet analysis agree well with the experimental results, particularly in the region with high frequencies.

As can be seen from figure 6, a zigzag mode appears in the experiment. Equation (4.22a) (zigzag mode) is very close to the experimental results, as shown in figure 7.

However, since the asymptotic expression is the first-order asymptotic result of the viscosity, the curve of the asymptotic analysis is slightly lower than that of the Floquet analysis. Furthermore, by adopting (4.17), we calculate the value of the second-order asymptotic and plot it in figure 7. Results suggest that at lower frequencies, the curves generated by the second-order asymptotic are more accurate than those generated by the first-order asymptotic because the second-order asymptotic includes more interfacial coupling effects. However, as shown in figure 7, at low frequencies corresponding to large wavenumbers, there are still discrepancies between theoretical and experimental results. Discrepancies between the results of these approaches decrease with increasing wavenumber (which is directly proportional to frequency). The observed difference between the experimental and theoretical critical vibration amplitudes in the low-frequency region can be attributed to the increase in importance of the viscous coupling effect as the frequency decreases.

6. Faraday instability modes

By analysing the relative displacements of the two interfaces, Pototsky & Bestehorn (2016) identify two types of modes: in-phase and anti-phase. The in-phase displacement corresponds to the zigzag (barotropic) mode, while the anti-phase displacement corresponds to the varicose mode. Nevertheless, the instability criterion of different modes is not well understood. In this section, we investigate further the Faraday instability modes from a theoretical point, and discuss the relationship between the interfacial coupling strength and the instability modes.

Equations (2.19) and (2.20) show that the coupling strength of two interfaces is determined by the dimensionless parameters $\beta _A$ and $\beta _B$, which range from 0 to 1. The region with strong coupling of two interfaces appears when $0 \ll \rho _A,\rho _B<1$, while the region with weak coupling of two interfaces happens when $\rho _A\sim 0$ and $\rho _B\sim 0$, corresponding to the case when (2.19), (2.20) are approximately decoupled.

6.1. Strongly coupled double-interface system

In this subsection, we perform Floquet analysis and find a direct solution of the Mathieu equations (2.19), (2.20) using the parameters at $f=60$ Hz from § 5. The initial conditions of the direct solution are set as $\xi _A = \xi _B = \xi _0$ and $\dot {\xi _A} = \dot {\xi _B} = 0$.

As shown in figure 8(a), the growth rate contours of the solution are calculated using the ODE89 function in MATLAB. The value of growth rate $s$ is obtained by fitting exponentially ($\exp (s_A t)$ and $\exp (s_B t)$) the envelope of the wave amplitude and then selecting the maximum value of the two interface growth rates for presentation. Furthermore, the Faraday tongue (black dotted line in figure 8a) is calculated by Floquet analysis, which is in agreement with the direct solution of the Mathieu equations, except for the narrow third tongue (pink dashed line in figure 8a) near the double-tongue stability curves. Similar to Pototsky & Bestehorn (2016), we also found two instability modes, the zigzag mode (point I), and the varicose mode (point II), as shown in figure 8(a). Two modes correspond to two extremes of double-tongued stability curves of the strongly coupled double-interface system. Figure 8(a) also shows that the zigzag mode is dominant in our experiments because of the lower Faraday threshold. Furthermore, we can use the zero-order approximate equation (4.4) to give an illustration of the occurrence of two instability modes. Substituting (4.7) into (4.4), we obtain

(6.1) \begin{equation} \left.\begin{gathered} \dfrac{\xi_A}{\xi_B} = \dfrac{\omega^2 \beta_A}{\omega^2-\varOmega_A^2},\\ \dfrac{\xi_A}{\xi_B} = \dfrac{\omega^2-\varOmega_B^2}{\omega^2 \beta_B}.\end{gathered}\right\} \end{equation}

Figure 8. Critical stability curve, growth rate and instability mode of the strongly coupled system and weakly coupled system: (a) strongly coupled double-interface system $h_2=1.262\ \mathrm {mm}$; (b) weakly coupled double-interface system $h_2=12.62\ \mathrm {mm}$; (c) typical modes and instabilities. The growth rate $s$ of waves is calculated by the ODE89 function in the strongly coupled double-interface system, the double-tongued stability curves of Floquet solutions, and the location of the three modes of instability. The black dotted line (the outline of the coloured area) is the first subharmonic stability curve calculated by Floquet analysis. The pink dashed line is the stability curve of the beating mode calculated by the ODE89 function.

The dispersion relation (4.11) can be obtained from (6.1). Furthermore, characteristic frequencies $\omega$ should satisfy $\omega >\max \{ \varOmega _A, \varOmega _B \}$ or $\omega <\min \{ \varOmega _A, \varOmega _B \}$ when $0<(\beta _A, \beta _B)<1$. We notice that when $\omega >\max \{ \varOmega _A, \varOmega _B \}$, $\xi _A/\xi _B>0$, that is, zigzag mode occurs. Similarly, when $\omega <\min \{ \varOmega _A, \varOmega _B \}$, $\xi _A/\xi _B<0$, namely, the varicose mode occurs. It should be noted that these two modes of instability do not alter the characteristics of subharmonic instability. As shown in figure 9(a), the Fourier analysis results show that the wave frequency is concentrated at $\omega$, the subharmonic frequency. Simultaneously, the triple frequency is observed, which can be explained by the first-order asymptotic approximation (4.12), (4.13). Although the results of the Floquet analysis are in good agreement with the direct solution of the Mathieu equations in most cases in figure 8(a), the ordinary differential equation solutions give a narrow third tongue near the double-tongued stability curves in figure 8(a). Figure 9(a) shows that the beating mode that exhibits primary frequencies is found at point III (red solid and dashed lines). Moreover, the average frequency of the beating mode is half of the external excitation frequency. Even though this phenomenon has not yet been found in experiments, we are still able to provide some results from a theoretical analysis. It can be seen from Appendix B that the two primary frequencies of the beating wave are the two characteristic frequencies, $\omega _{0,1}$ and $\omega _{0,2}$, obtained from the dispersion relation. The centre frequency of the beating wave is $(\omega _{0,1}+\omega _{0,2})/2$ in figure 8(a). Figure 9(a) confirms that the frequency difference is $| \omega _{0,1}-\omega _{0,2} |$, and this phenomenon also exists near the triple frequency. The beating mode obtained here is similar to the combination resonances discovered by Kidambi (2013). Both are generated by the coupling effects of the system. However, the beating instability involves a dynamic coupling between two interfaces, whereas the coupling involves interactions between different modes induced by a fixed contact line.

Figure 9. Results of fast Fourier transform in strongly coupled and weakly coupled double-interface systems: (a) strongly coupled double-interface system, $h_2=1.262\ \mathrm {mm}$; (b) weakly coupled double-interface system, $h_2=12.62\ \mathrm {mm}$. The solid line refers to interface $A$, while the dotted line refers to interface $B$, and $\varOmega _w$ is the frequency of the wave.

6.2. Weakly coupled double-interface system

To discuss the weakly coupled double interface system, we changed the thickness of the sodium hydroxide solution to 12.62 mm, while keeping the thickness of the GaInSn and the properties of the liquid constant. It has been found by Ward et al. (2019b) that if the middle fluid layer is thick enough, then the stability curve of the three-fluid system is composed primarily of two separate two-fluid systems. Equations (4.4) are approximately decoupled when the coupling between the two interfaces is weak ($\beta _A,\beta _B \sim 0$). Meanwhile, the dispersion relation (4.11) is simplified as (4.23a,b).

According to figure 8(b), $\varOmega _A$ and $\varOmega _B$ are the characteristic frequencies of interface $A$ and interface $B$, respectively. Here, $\varOmega _A$ is determined by the GaInSn/solution interface, while $\varOmega _B$ is determined by the solution/air interface. In other words, $\omega _{0,A} = \varOmega _A$ and $\omega _{0,B} = \varOmega _B$, where $\omega _{0,A}$ and $\omega _{0,B}$ are the frequencies of the $A$-interface instability (point IV) and $B$-interface instability (point V), respectively. In addition, as shown in figure 8(b), the critical stability curve can be separated into the curves of interfaces $A$ and $B$, due to (2.17). Therefore, at the extreme points of these two instability curves, the instability modes correspond to the instabilities of interfaces $A$ and $B$, respectively. At a codimension-two point in figure 8(b), the instability of the two interfaces ($A$&$B$ instability at point VI) with different wavenumbers appears simultaneously, which agrees well with that from Ward et al. (2019b). Figure 9(b) shows that the frequency of waves is concentrated at $\omega$ in the weakly coupled double-interface system.

6.3. Relationship between coupling strength and instability modes

The previous results indicate that changes in coupling strength affect the Faraday instability mode. Under strongly coupled conditions ($0 \ll \beta _A,\beta _B<1$), zigzag mode, varicose mode, and beating mode occur (§ 6.1), while under weakly coupled conditions ($\beta _A,\beta _B \sim 0$), the instabilities of the two interfaces are independent of each other. All instability modes, except the beating mode, satisfy approximately the dispersion relation (4.11). Equation (4.11) gives us two solutions, $\omega _{0,1}$ and $\omega _{0,2}$:

(6.2)\begin{equation} \left.\begin{gathered} \omega_{0,1}=\sqrt{\frac{\varOmega_A^2+\varOmega_B^2+\sqrt{\left(\varOmega_A^2- \varOmega_B^2\right)^2+4 \beta_A \beta_B \varOmega_A^2 \varOmega_B^2}}{2\left(1- \beta_A \beta_B\right)}},\\ \omega_{0,2}=\sqrt{\frac{\varOmega_A^2+\varOmega_B^2-\sqrt{\left(\varOmega_A^2-\varOmega_B^2\right)^2+4 \beta_A \beta_B \varOmega_A^2 \varOmega_B^2}}{2\left(1- \beta_A \beta_B\right)}}. \end{gathered}\right\} \end{equation}

In figure 10(b), data of mode $(3, 4)$ are obtained by substituting (6.2) into any equation in (6.1). As $kh_2$ increases, the zigzag mode transforms gradually into $B$-interface instability, while the varicose mode transitions gradually to $A$-interface instability. This change is caused primarily by the coupling strength of the two interfaces. It can be seen that different instability modes are specific manifestations of the two solutions of the dispersion relation under different coupling degrees. To describe the strength of interface coupling, the parameters $\beta _A$ and $\beta _B$, which come from the derivation of the equation, are used. It can be seen from figure 10(a) that $\beta _A,\beta _B$ decreases rapidly with the increase of $kh_2$. The influences of $\beta _A$ and $\beta _B$ are also reflected in the dispersion relationships, as shown in figure 10(c). As $kh_2$ increases, $\omega _{0,1}$ approaches $\varOmega _A$, and $\omega _{0,1}$ approaches $\varOmega _B$.

Figure 10. (a) The relationship between coupling strength coefficient $\beta _A, \beta _B$ and dimensionless number $kh_2$ of mode $(3,4)$, where $h_1=12.63\ \textrm {mm}$, $R=27.5\ \textrm {mm}$, $\rho _1 = 6360\ \textrm {kg}\ \textrm {m}^{-3}$ and $\rho _2 = 1041\ \textrm {kg}\ \textrm {m}^{-3}$, calculated by (2.14a,b). (b) The relationship between the amplitude ratio of the two interfaces $\xi _A / \xi _B$ and dimensionless number $kh_2$ of mode $(3, 4)$, calculated by (6.1) and (6.2). (c) The two solid lines are calculated from (6.2), and the two dashed lines are calculated from (2.14a,b), for mode $(3,4)$.

However, in figure 11(b), as $kh_2$ increases for mode $(4, 4)$, zigzag mode transforms gradually to $A$-interface instability, and varicose mode transits to $B$-interface instability. Moreover, figure 11(c) shows that $\omega _{0,1}$ converges to $\varOmega _A$, and $\omega _{0,2}$ converges to $\varOmega _B$ when $\varOmega _A > \varOmega _B$. It is observed in figures 10(c) and 11(c) that the modal conversion directions between mode $(3, 4)$ and mode $(4, 4)$ are related to $\varOmega _A/\varOmega _B$. Assuming $h_j \to \infty$, $k_c=\sqrt {2g( \rho _2^2-\rho _1 \rho _3 )/(\sigma _A( \rho _2+\rho _3 )-\sigma _B( \rho _1+\rho _2 ))}$ is obtained from (2.14a,b) when $\varOmega _A/\varOmega _B=1$, which is the critical wavenumber of the modes transition. From figure 12, when $k< k_c$, $\omega _{0,1}$ converges to $\varOmega _B$, and $\omega _{0,2}$ converges to $\varOmega _A$, similar to the previous discussion of mode $(3,4)$. When $k>k_c$, $\omega _{0,1}$ converges to $\varOmega _A$, and $\omega _{0,2}$ converges to $\varOmega _B$, similar to the previous discussion of mode $(4,4)$.

Figure 11. (a) The relationship between coupling strength coefficient $\beta _A, \beta _B$ and dimensionless number $kh_2$ of mode $(4,4)$, where $h_1=12.63\ \textrm {mm}$, $R=27.5\ \textrm {mm}$, $\rho _1 = 6360\ \textrm {kg}\ \textrm {m}^{-3}$ and $\rho _2 = 1041\ \textrm {kg}\ \textrm {m}^{-3}$, calculated by (2.14a,b). (b) The relationship between the amplitude ratio of the two interfaces $\xi _A / \xi _B$ and dimensionless number $kh_2$ of mode $(4,4)$, calculated by (6.1) and (6.2). (c) The two solid lines are calculated from (6.2), and the two dashed lines are calculated from (2.14a,b) for mode $(4,4)$.

Figure 12. The curves of $\varOmega _A/\varOmega _B$, $\omega _{0,1}/\varOmega _A$, $\omega _{0,1}/\varOmega _B$, $\omega _{0,2}/\varOmega _A$ and $\omega _{0,2}/\varOmega _B$ converge at the point $(k_c,1)$, where $R=27.5\ \textrm {mm}$, $\rho _1 = 6360\ \textrm {kg}\ \textrm {m}^{-3}$ and $\rho _2 = 1041\ \textrm {kg}\ \textrm {m}^{-3}$, calculated by (2.14a,b). The wavenumber of mode $(3,4)$ is less than $k_c$, and the wavenumber of mode $(4,4)$ is greater than $k_c$.

7. Conclusion

By applying the viscous model derived by Labrador et al. (2021), the present paper examines Faraday instability of two-interface fluid systems from both theoretical and experimental perspectives to study the effects caused by the coupling of two interfaces. The Mathieu equations are expanded asymptotically. The dispersion relation is obtained through the zeroth-order approximation expansion, while the critical vibration amplitude for the onset of Faraday instability is obtained through the second-order approximation expansion. These results are then compared to the results from the Floquet analysis of Mathieu equations and the Faraday experimental system. The dispersion relation shows good agreement between experimental data and theoretical predictions for the zigzag mode. The results of the Floquet analysis are in good agreement with the experiments regarding the critical vibration amplitude for the onset of the Faraday instability, but there are differences between the theoretical and experimental results in the low-frequency region corresponding to the large wavenumber. Furthermore, both theoretical and experimental results demonstrate that the discrepancy between the outcomes of the two approaches diminishes as the wavenumber increases (wavenumber is directly proportional to frequency). The critical oscillation amplitude in the low-frequency region differs from the theoretical analysis due to the damping effects caused by wall and interface viscous shear (Batson et al. 2013). However, the difference between Floquet and asymptotic analyses lies in their respective accuracies: the asymptotic analysis is a second-order approximation, while the Floquet series is truncated at the 40th order.

Furthermore, several Faraday instability modes are discussed. According to the Mathieu equations, fluid systems can be classified as strongly coupled or weakly coupled. In strongly coupled systems, the zigzag and varicose modes can be obtained directly by solving the ordinary differential equation. These two types of instability correspond to the two extreme points of the instability tongue. The instability criterion for these two modes can be determined further by using the zeroth-order approximate expansion of the Mathieu equations, $\omega _{{zigzag}}>\max \{ \varOmega _A,\varOmega _B \}$ and $\omega _{{varicose}}<\min \{ \varOmega _A,\varOmega _B \}$. In addition, at the extreme point of the third tongue (point III in figure 8a), we can find the beating mode and obtain the centre frequency $(\omega _{0,1}+\omega _{0,2})/2$ and frequency difference $|\omega _{0,1}-\omega _{0,2}|$. In weakly coupled fluid systems, the motions of interfaces $A$ and $B$ are decoupled. At the extreme point of the critical instability curve near $\omega =\varOmega _A$, interface $A$ becomes unstable. Similarly, interface $B$ becomes unstable at the extreme point of the critical instability curve near $\omega =\varOmega _B$. At a codimension-two point, where two critical instability curves of interfaces $A$ and $B$ intersect, $A$&$B$ instability occurs. Furthermore, under different coupling conditions, the instability modes are special manifestations of two solutions of the dispersion equation. For standing wave modes with wavenumbers less than the critical wavenumber $k_c$, the zigzag mode transforms gradually into the $B$-interface instability mode, and the varicose mode transforms gradually into the $A$-interface instability mode, as $\beta _A$ and $\beta _B$ increase. For standing wave modes with wavenumbers greater than $k_c$, the zigzag mode transforms gradually into the $A$-interface instability mode, and the varicose mode transforms gradually into the $B$-interface instability mode.

The current study provides valuable insights for further investigations of the dynamics of multi-layer fluid interfaces. By deriving the dispersion relation of the three-layer fluid system, we provide a theoretical reference for analysing interface instability problems in various systems, including liquid metal batteries (Horstmann, Weber & Weier 2018). Moreover, we propose a theoretical framework for understanding multi-layer fluid systems with low viscosity. Nevertheless, further research is necessary to address the interface coupling issue in high-viscosity fluid systems.