2. Experimental platform

Circular orbit oscillations induce interface waves in a viscous fluid, allowing us to investigate the damping mechanisms of multi-layer stratified fluid interface waves. Similar to the experimental set-up of Horstmann et al. (Reference Horstmann, Wylega and Weier2019), as shown in figure 1, we construct a physical test system with three-layer fluid. A cylindrical container of inner diameter $2R = 80$ mm and height $H = 100$ mm is filled to the brim with three immiscible liquids in layers. The bottom layer (fluid 1) thickness is fixed at $h_1$ = 30 mm with the thickness $h_2$ of the middle layer (fluid 2) and the thickness $h_3 = H-h_1-h_2$ of the top layer (fluid 3). The three-layer fluid system (denoted as $j=1, 2, 3$ for fluids 1, 2 and 3, respectively), with densities $\rho _j$ , kinematic viscosities $\nu _j$ and thicknesses $h_j$ , is supported by an aluminium frame and fixed on a Kuhner ES-X lab-shaker (500 mm $\times$ 420 mm). The diameter of the shaker, ${d}_s$ , is 25 mm. The angular frequency $f = \varOmega /2\pi$ . In the cylindrical coordinate system $O(r, \varphi , z)$ , gravity is represented by the acceleration, $g$ , acting in the negative $z$ -direction. The interface tension coefficient of interface $A$ between fluids 1 and 2 is $\sigma _A$ at $z=\xi _A ( r,\varphi ,z )$ . Similarly, the interface tension coefficient of interface $B$ between fluids 2 and 3 is $\sigma _B$ at $z=\xi _B ( r,\varphi ,z )$ . Following the assumption suggested by Horstmann et al. (Reference Horstmann, Herreman and Weier2020), we neglect capillary forces acting on the contact line and assume a 90-degree contact angle between the interface and the cylindrical wall. The origin of the moving coordinate system $O(r, \varphi , z)$ is defined on a circle with a fixed point ${O}'$ as the centre and ${d}_s$ as the diameter. In the ${O}'({r}' , {\varphi }' ,{z}')$ coordinate system, point $O$ moves uniformly around point ${O}'$ . The velocity vector $V_0$ describing the orbital motion of the container can be expressed in the cylindrical coordinate system of the inertial reference frame ${O}'(r', {\varphi }', z')$ as $\boldsymbol{V}_0 ( t ) = ( -{d}_s\varOmega \sin ( \varOmega t-{\varphi }' )/2, {d}_s\varOmega \cos ( \varOmega t-{\varphi }' )/2,0 )$ .

An Air∣P∣S system (i.e. of air/paraffin oil/silicone oil) is selected to measure simultaneously the amplitude and phase difference of the two fluid interfaces. The interfacial tension coefficient between air and paraffin oil is approximately 20 mN m^–1. The interfacial tension coefficient between paraffin oil and silicone oil is approximately 5 mN m^–1. The ultrasonic technique to measure the positions of the interfaces is based on time to flight. Further, we extend this measurement method to two fluid interfaces (Appendix A). Figure 1(a) shows that an ultrasonic probe with a centre frequency of 4 MHz and a diameter of 8 mm is positioned on the bottom surface at $r_0$ = 30 mm. We use multichannel pulsed ultrasonic Doppler velocimetry (MPUDV) (Pan et al. Reference Pan, Huang, Lyu, Yang and Ni2025) supplied by the University of Chinese Academy of Science (UCAS) for the acoustic measurements of interface wave motion. As shown in figure 1(a), the elevation of interface $A$ (the interface between fluid 1 and fluid 2) above the base of the container is $\xi _A\approx c_1 \Delta t_A/2$ , where $c_1$ represents the speed of sound in fluid 1 and $\Delta t_A$ is the time to travel from transmission to reception at interface $A$ . The elevation of interface $B$ above the base of the container is $\xi _B\approx \xi _A + c_2 (\Delta t_B-\Delta t_A )/2$ , where $c_2$ is the speed of sound in fluid 2 and $\Delta t_B$ is the time to travel from transmission to reception at interface $B$ . The amplitudes, $\delta _A$ and $\delta _B$ , and the phase difference, $\Delta \varphi$ , of interface waves are obtained by ultrasonic measurement.

In addition, we select a P∣S∣G system (i.e. of paraffin oil/silicone oil/glycerol) to use the BOS method for measurements in our experiments. The interfacial tension coefficient between silicone oil (S) and glycerol (G) is approximately 18.2 mN m^–1. BOS is an optical method designed to measure the instantaneous topography of the interface between two transparent fluids, which has been applied to measure interface waves. The visual displacement of random speckle patterns through interface refraction is the theoretical basis of this method. The surface gradient field is determined by the apparent displacement between the refracted image and a reference image taken when the surface is flat. Following the work of Moisy, Rabaud & Salsac (Reference Moisy, Rabaud and Salsac2009), we can show that the virtual displacement $\delta \boldsymbol{r}$ reflects a weighted sum of the gradients at the two interfaces (Appendix B), $ ( 1-n_2/n_1 )h_1 \triangledown \xi _A + ( 1-n_3/n_2 ) ( h_1 +h_2 ) \triangledown \xi _B \sim -\delta \boldsymbol{r}$ , where $n_j$ represents the refractive index of the fluid $j$ . As shown in figure 1, a diffuse backlight LED array ( $150\,{\textrm{mm}}\times 150\,{\textrm{mm}}$ ) is placed below the container. Additionally, a random dot pattern of approximately 20 dots mm^– $^2$ , generated using the PIVMat Toolbox, is placed between the container and the LED array. We captured images of the pattern at a resolution of 0.06 mm pixel^–1 with an industrial camera (U3P630-H), positioned approximately 200 mm above the container at a frame rate of 70 fps. The system is placed on the shaker such that the camera and light source moved synchronously with the experimental set-up. Finally, digital image correlation (DIC) is performed on the point patterns captured during motion and the stationary point patterns using MATLAB.

3. Theory

The incompressible N-S equations and boundary conditions for the viscous interface wave mechanics of three-layer fluid systems are listed in Appendix C in terms of velocity $\boldsymbol{v}_j$ , pressure $p_j$ , density $\rho _j$ and kinematic viscosity $\nu _j$ $(j=1,2,3)$ . We define a small parameter $\epsilon = 2Fr$ , where $Fr = {d}_s \varOmega ^2/g$ is the Froude number, and perform a perturbation-based asymptotic expansion as follows:

(3.1) \begin{equation} \begin{aligned} \boldsymbol{v}_j &=\boldsymbol{u}_j+{\textrm{O}}\big(\epsilon ^2\big ), \\p_j &= p_{j, 0}+\rho _j \frac {{d}_s \varOmega ^2 r}{2} \cos (\varOmega t-\varphi )+\hat {p}_j+{\textrm{O}}\big (\epsilon ^2\big ), \\\xi _A &= h_1+\eta _A+{\textrm{O}}\big (\epsilon ^2\big ), \\\xi _B &= h_1+h_2+\eta _B+{\textrm{O}}\big(\epsilon ^2\big). \end{aligned} \end{equation}

From the perturbation expansion, we obtain the first-order incompressible N-S equations (3.2b ),

(3.2a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_j &=0, \end{align}

(3.2b) \begin{align} \frac {\partial \boldsymbol{u}_j}{\partial t} &=-\frac {1}{\rho _j} \boldsymbol{\nabla }\hat {p}_j+v_j {\nabla}^2 \boldsymbol{u}_j, \end{align}

and corresponding boundary conditions (3.3),

(3.3a) \begin{align} \textrm{Wall}: &\; \left .\boldsymbol{u}_j\right |_{z=0\|\| z=\sum _1^3 h_j \| r=R}=0, \end{align}

(3.3b) \begin{align} \textrm{Interface}\; A: \;& \left .\boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1}=\left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1}=\frac {\partial \eta _A}{\partial t}, \end{align}

(3.3c) \begin{align} & \left .\boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1}=\left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1}\!, \end{align}

(3.3d) \begin{align} & \left .\boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1}=\left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1}, \end{align}

(3.3e) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{1, l s}-\tau _{2, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1}=\sigma _A {\nabla}^2 \eta _A, \end{align}

(3.3f) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{1, l s}-\tau _{2, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1}=0, \end{align}

(3.3g) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{1, l s}-\tau _{2, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1}=0, \end{align}

(3.3h) \begin{align} \textrm{Interface}\; B: &\; \left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1+h_2}=\left .\boldsymbol{u}_3 \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1+h_2}=\frac {\partial \eta _B}{\partial t}, \end{align}

(3.3i) \begin{align} & \left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1+h_2}=\left .\boldsymbol{u}_3 \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1+h_2}\!, \end{align}

(3.3j) \begin{align} & \left .\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1+h_2}=\left .\boldsymbol{u}_3 \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1+h_2}\!, \end{align}

(3.3k) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{2, l s}-\tau _{3, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_z\right |_{z=h_1+h_2}=\sigma _B {\nabla}^2 \eta _B, \end{align}

(3.3l) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{2, l s}-\tau _{3, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_x\right |_{z=h_1+h_2}=0, \end{align}

(3.3m) \begin{align} & \left .\boldsymbol{e}_z \boldsymbol{\cdot }\left (\tau _{2, l s}-\tau _{3, l s}\right ) \boldsymbol{\cdot }\boldsymbol{e}_y\right |_{z=h_1+h_2}=0, \end{align}

with $\tau _{j,ls}$ as the stress tensor of fluid $j$ .

Table 1.

Experimental parameters of the two sets of experiments.

3.1. The modified model

In practice, it is not easy to find analytical solutions for the dynamics of viscous waves. In the past, potential theory was applied in many studies to approximate the velocity field based on inviscid, irrotational and incompressible flow assumptions. Nevertheless, potential theory has successfully predicted non-resonant sloshing (Ibrahim Reference Ibrahim2005). As described by Faltinsen & Timokha (Reference Faltinsen and Timokha2009), although potential flow theory is strictly only applicable to irrotational flow, it allows the inclusion of viscous damping rates in the model equations. In this way, provided the wave amplitude is sufficiently small, the modified model accounts for linear and weakly nonlinear phenomena (Faltinsen & Timokha Reference Faltinsen and Timokha2009). Horstmann et al. (Reference Horstmann, Herreman and Weier2020) proposed a corresponding modified model in a two-layer viscous fluid system, adding the explicit viscous damping rates in cylindrical containers developed by Herreman et al. (Reference Herreman, Nore, Guermond, Cappanera, Weber and Horstmann2019). Similarly, a corresponding modified model can also be derived based on explicit damping rates for three-layer viscous fluid cylindrical containers, as presented by Herreman et al. (Reference Herreman, Wierzchalek, Horstmann, Cappanera and Nore2023). However, the coupling effects between two fluid interfaces need to be evaluated for a three-layer fluid system. We established this model based on Horstmann’s theoretical work (Horstmann et al. Reference Horstmann, Herreman and Weier2020, Reference Horstmann, Anders, Kelley and Weier2021; Horstmann Reference Horstmann2021).

3.1.1. Potential theory

Based on potential flow theory (ignoring viscosity), the velocity field $\boldsymbol{u}_j$ can be represented by a potential function $\varPhi _j$ . The velocity potential and interface wave functions should have the following functional forms to satisfy the slip condition at walls and interfaces:

(3.4a) \begin{align} \eta _A(r, \phi ) & =\sum _{n=1}^{\infty } \beta _{A, n}(\phi ) \frac {J_1\big (k_{1, n} r\big )}{J_1\big (k_{1, n} R\big )}, \end{align}

(3.4b) \begin{align} \eta _B(r, \phi ) & =\sum _{n=1}^{\infty } \beta _{B, n}(\phi ) \frac {J_1\big (k_{1, n} r\big )}{J_1\big (k_{1, n} R\big )}, \end{align}

(3.4c) \begin{align} \varPhi _1(r, \phi , z) & =\sum _{n=1}^{\infty } \frac {\varOmega \dot {\beta }_{A, n}}{k_{1, n}} \frac {\cosh k_{1, n} z}{\sinh k_{1, n} h_1} \frac {J_1\big (k_{1, n} r\big )}{J_1\big (k_{1, n} R\big )}, \end{align}

(3.4d) \begin{align} &{\varPhi _2(r, \phi , z) = -\! \sum _{n=1}^{\infty} \! \left ( \! \frac {\varOmega \dot {\beta }_{A, n}}{k_{1, n}} \frac {\cosh k_{1, n} (z \!- \!h_1 \!- \!h_2)}{\sinh k_{1, n} h_2} \!-\! \frac {\dot {\varOmega \beta }_{B, n}}{k_{1, n}} \frac {\cosh k_{1, n} (z \!- \!h_1)}{\sinh k_{1, n} h_2}\! \right ) \!\!\!\! \quad \frac {J_1 (k_{1, n} r)}{J_1 (k_{1, n} R)}, }\end{align}

(3.4e) \begin{align} \varPhi _3(r, \phi , z) & =-\sum _{n=1}^{\infty } \frac {\varOmega \dot {\beta }_{B, n}}{k_{1, n}} \frac {\cosh k_{1, n}\left (z-h_1-h_2-h_3\right )}{\sinh k_{1, n} h_3} \frac {J_1\big (k_{1, n} r\big )}{J_1\big (k_{1, n} R\big )}, \end{align}

with, $k_{1,n}=\varepsilon _{1,n} / R, J_1^{\prime } (\varepsilon _{1,n} )=0, \phi =\varOmega t-\varphi$ , where the quantities $\eta _A$ and $\eta _B$ denote the first-order approximations of the interfacial waves at interfaces A and B, respectively, with $J_1(k_{1,n} r)$ representing the Bessel function of the first kind of order one. The interfacial dynamic boundary conditions (3.3e , 3.3k ) may be simplified to give

(3.5a) \begin{align} & \left .\rho _1 \varOmega \frac {\textrm{d} \varPhi _1}{\textrm{d} \phi } \right |_{z=h_1}-\left .\rho _2 \varOmega \frac {\textrm{d} \varPhi _2}{\textrm{d} \phi }\right |_{z=h_1} +\big ((\rho _1-\rho _2) g-\sigma _A {\nabla} _{(r, \varphi )}^2 \big ) \nonumber\\ & \eta _A = (\rho _1-\rho _2) \frac {{d}_s \varOmega ^2 r}{2} \cos \phi , \end{align}

(3.5b) \begin{align} & \left .\rho _2 \varOmega \frac {\textrm{d} \varPhi _2}{\textrm{d} \phi }\right |_{z=h_1+h_2}-\left .\rho _3 \varOmega \frac {\textrm{d} \varPhi _2}{\textrm{d} \phi }\right |_{z=h_1+h_2} +\big ((\rho _2-\rho _3) g-\sigma _B {\nabla} _{(r, \varphi )}^2\big)\nonumber\\ & \eta _B =\left (\rho _2-\rho _3\right ) \frac {{d}_s \varOmega ^2 r}{2} \cos \phi , \end{align}

with ${\nabla} _{(r, \varphi )}^2 =( {\partial ^2 }/{\partial r^2})+ ({1}/{r}) ({\partial }/{\partial r})+ ({1}/{r^2}) ({\partial ^2 }/{\partial \theta ^2})$ . Substituting (3.4) into (3.5), we obtain

(3.6a) \begin{align} & \frac {J_1\left (k_{1, n} r\right )}{J_1\left (k_{1, n} R\right )}\left (\varOmega ^2 \frac {\textrm{d}^2 \beta _{A, n}}{\textrm{d} \phi ^2} -\rho _{A, n} \varOmega ^2 \frac {\textrm{d}^2 \beta _{B, n}}{\textrm{d} \phi ^2}+\varOmega _{A, n}^2 \beta _{A, n}\right )\nonumber\\ &\quad =\frac {{d}_s}{2} \frac {\left (\rho _1-\rho _2\right ) k_{1 n} \varOmega ^2 r}{\rho _1 \operatorname {coth} (k_{1 n} h_1)+\rho _2 \operatorname {coth} (k_{1 n} h_2)} \sin \phi , \end{align}

(3.6b) \begin{align} & \frac {J_1\left (k_{1, n} r\right )}{J_1\left (k_{1, n} R\right )}\left (\varOmega ^2 \frac {\textrm{d}^2 \beta _{B, n}}{\textrm{d} \phi ^2}-\rho _{B, n} \varOmega ^2 \frac {\textrm{d}^2 \beta _{A, n}}{\textrm{d} \phi ^2}+\varOmega _{B, n}^2 \beta _{B, n}\right ) \nonumber\\ &\quad =\frac {{d}_s}{2} \frac {\left (\rho _2-\rho _3\right ) k_{1 n} \varOmega ^2 r}{\rho _2 \operatorname {coth} (k_{1 n} h_2)+\rho _3 \operatorname {coth} (k_{1 n} h_3)} \sin \phi , \end{align}

with

(3.7) \begin{eqnarray} \rho _{A, n}=\frac {\rho _2 \operatorname {csch} (k_{1,n} h_2)}{\rho _1 \operatorname {coth} (k_{1,n} h_1)+\rho _2 \operatorname {coth} (k_{1,n} h_2)}, \rho _{B, n} = \frac {\rho _2 \operatorname {csch} (k_{1,n} h_2)}{\rho _2 \operatorname {coth} (k_{1,n} h_2)+\rho _3 \operatorname {coth} (k_{1,n} h_3)}, \nonumber\\[8pt] \varOmega _{A, n}^2 = \frac {(\rho _1-\rho _2)g k_{1,n}+\sigma _A k_{1,n}^3}{\rho _1 \operatorname {coth} (k_{1,n} h_1)+\rho _2 \operatorname {coth} (k_{1,n} h_2)}, \varOmega _{B, n}^2 = \frac {(\rho _2-\rho _3)g k_{1,n}+\sigma _B k_{1,n}^3}{\rho _2 \operatorname {coth} (k_{1,n} h_2)+\rho _3 \operatorname {coth} (k_{1,n} h_3)}, \!\!\!\!\nonumber\\\end{eqnarray}

where $\varOmega _{A,n}$ and $\varOmega _{B,n}$ represent the natural frequencies (Huang et al. Reference Huang, Zhuo, Yang and Ni2024) associated with interface $A$ and interface $B$ , respectively. The parameters $\rho _{A,n}$ and $\rho _{B,n}$ quantify the coupling strength between the two interfacial waves. When the shallow assumption $k_{1,n} h_j\ll 1$ is satisfied,

(3.8) \begin{equation} \rho _{A, n} \approx \frac {1}{1+\left (h_2 / h_1\right )\left (\rho _1 / \rho _2\right )},\quad \rho _{B, n} \approx \frac {1}{1+\left (h_2 / h_3\right )\left (\rho _3 / \rho _2\right )} \end{equation}

and

(3.9) \begin{equation} \frac {\rho _{A, n}}{\rho _{B, n}} \approx \frac {1+\left (\rho _3 / \rho _2\right )\left (h_2 / h_3\right )}{1+\left (\rho _1 / \rho _2\right )\left (h_2 / h_1\right )}. \end{equation}

Figure 2.

Variation of the coupling coefficient $\rho _{A, n} / \rho _{B, n}$ with $k_{1, n} h_2$ in the Air∣P∣S system.

Furthermore, when $h_2 \rightarrow 0$ , $\lim _{h_2 \rightarrow 0} ({\rho _{A, n}}/{\rho _{B, n}}) \approx 1$ is obtained, indicating a strong interfacial coupling effect.

When the deep assumption $k_{1, n} h_j \gg 1$ is satisfied,

(3.10) \begin{equation} \rho _{A, n} \approx 0, \,\, \rho _{B, n} \approx 0 \end{equation}

and

(3.11) \begin{equation} \frac {\rho _{A, n}}{\rho _{B, n}} \approx \frac {1+\rho _3 / \rho _2}{1+\rho _1 / \rho _2}, \end{equation}

which indicate a weak interfacial coupling effect.

Therefore, the closer $\rho _{A, n} / \rho _{B, n}$ is to one, the stronger the interfacial coupling effect; conversely, as $\rho _{A, n} / \rho _{B, n}$ approaches $ (1+\rho _3 / \rho _2 ) / (1+\rho _1 / \rho _2 )$ , the coupling becomes weaker. To illustrate this intuitively, we take the Air∣P∣S system as an example and plot the variation of $\rho _{A, n} / \rho _{B, n}$ concerning $k_{1, n} h_2$ in figure 2.

From the orthogonality of the Bessel function $J_1 (k_{1,n} r )$ , (3.6) is reduced to

(3.12a) \begin{align} \frac {\textrm{d}^2 \beta _{A, n}}{\textrm{d} \phi ^2} -\rho _{A, n} \frac {\textrm{d}^2 \beta _{B, n}}{\textrm{d} \phi ^2} +\frac {\varOmega _{A, n}^2}{\varOmega ^2} \beta _{A, n}&=F_{A,n} {d}_s \sin \phi , \end{align}

(3.12b) \begin{align} \frac {\textrm{d}^2 \beta _{B, n}}{\textrm{d} \phi ^2}-\rho _{B, n} \frac {\textrm{d}^2 \beta _{A, n}}{\textrm{d} \phi ^2}+\frac {\varOmega _{B, n}^2}{\varOmega ^2} \beta _{B, n}&=F_{B,n} {d}_s \sin \phi , \end{align}

where

(3.13) \begin{equation} \begin{aligned} F_{A, n}=\frac {\left (\rho _1-\rho _2\right ) k_{1,n} R}{\left (\rho _1 \operatorname {coth} k_{1,n} h_1+\rho _2 \operatorname {coth} k_{1,n} h_2\right )\big (\varepsilon _{1,n}^2-1\big )J_1\left (k_{1, n} R\right )}, \\[6pt] F_{B, n}=\frac {\left (\rho _2-\rho _3\right ) k_{1,n} R}{\left (\rho _2 \operatorname {coth} k_{1,n} h_2+\rho _3 \operatorname {coth} k_{1,n} h_3\right )\big (\varepsilon _{1,n}^2-1\big )J_1\left (k_{1, n} R\right )}. \end{aligned} \end{equation}

The coefficients $F_{A,n}$ and $F_{B,n}$ quantify the influence of density stratification on the wave amplitudes at interfaces A and B, respectively, in a cylindrical container with finite fluid depth. As $h_j \to \infty$ , $F_{A,n}$ and $F_{B,n}$ can be approximated as

(3.14) \begin{equation} \begin{aligned} F_{A, n} \approx \frac {A_{1,2} k_{1, n} R}{\big (\varepsilon _{1, n}^2-1\big ) J_1\big (k_{1, n} R\big )}, \quad F_{B, n} \approx \frac {A_{2,3} k_{1, n} R}{\big (\varepsilon _{1, n}^2-1\big ) J_1\big (k_{1, n} R\big )}, \end{aligned} \end{equation}

where $A_{1,2}=(\rho _1-\rho _2)/(\rho _1+\rho _2)$ denotes the Atwood number between fluid 1 and fluid 2, and $A_{2,3}=(\rho _2-\rho _3)/(\rho _2+\rho _3)$ denotes the Atwood number between fluid 2 and fluid 3.

3.1.2. Viscous damping

Herreman et al. (Reference Herreman, Wierzchalek, Horstmann, Cappanera and Nore2023) carried out a detailed study of the viscous damping terms in a three-layer fluid system, starting with the energy equation,

(3.15) \begin{equation} \begin{aligned} &\frac {\textrm{d}}{\textrm{d} t} \underbrace {\sum _{j=1,2,3} \frac {1}{2} \int _{V_j} \rho _j\left \|\boldsymbol{u}_j\right \|^2 \, \textrm{d} V}_{E_c}+\\ &\frac {\textrm{d}}{\textrm{d} t} \underbrace {\int _{S_A}\left [\frac {1}{2}\left (\rho _2-\rho _1\right ) g \xi _A^2+\sigma _A \sqrt {1+\left \|\boldsymbol{\nabla }\xi _A\right \|^2}\right ] \textrm{d} S+\int _{S_B}\left [\frac {1}{2}\left (\rho _3-\rho _2\right ) g \xi _B^2+\sigma _B \sqrt {1+\left \|\boldsymbol{\nabla }\xi _B\right \|^2}\right ] \textrm{d} S}_{E_p}=\\ & \qquad\underbrace {-\sum _{j=1,2,3} 2 \rho _j v_j \int _{V_j} {\varepsilon }_j: {\varepsilon }_j \,\textrm{d} V}_D , \end{aligned} \end{equation}

with kinetic energy $E_c$ , potential energy $E_p$ , viscous dissipation energy $D$ and the shear rate tensor $\boldsymbol{\varepsilon }_j= (\partial u_{j, m} / \partial x_n+\partial u_{j, n} / \partial x_m ) / 2$ . The simplification of viscous dissipation energy $D$ leads to the following expression:

(3.16) \begin{equation} \begin{aligned} D=\underbrace {-\sum _{j=1,2,3} \rho _j v_j \int _{V_j}\left \|\boldsymbol{\nabla }\times \boldsymbol{u}_j\right \|^2 \, \textrm{d} V}_{D^{\prime }} \underbrace {-2 \sum _{j=1,2,3} \rho _j v_j \int _{V_j} \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\left (\boldsymbol{u}_j \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u}_j\right )\, \textrm{d} V}_{D^{\prime \prime }}. \end{aligned} \end{equation}

In this context, the main contributions come from the boundary layer ( $D^{\prime }$ ) and the fluid bulk ( $D^{\prime \prime }$ ). The viscous attenuation surface wave field and the interface displacement are $\boldsymbol{u}_j \sim \textrm{e}^{-\lambda t}, \eta \sim \textrm{e}^{-\lambda t}$ .

According to Herreman et al. (Reference Herreman, Wierzchalek, Horstmann, Cappanera and Nore2023), it can be observed that viscous damping consists of three components: wall damping $\lambda ^{{BL}}$ , interface damping $\lambda ^{{IL}}$ and internal fluid damping $\lambda ^{{Int}}$ . Here, $\lambda =\lambda ^{{BL}}+\lambda ^{{IL}}+\lambda ^{{Int}}$ , with $\lambda ^{{BL}},\lambda ^{{IL}}\sim k$ and $\lambda ^{{Int}}\sim k^2$ . Kumar & Tuckerman (Reference Kumar and Tuckerman1994) have provided a solution for $\lambda ^{{Int}}$ . As for the viscous dissipation energy in the fluid body, $D^{\prime \prime }$ can be approximated as

(3.17) \begin{equation} \begin{aligned} D^{\prime \prime }&=-2 \sum _{j=1,2,3} \rho _j v_j \int _{V_j} \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\left (\boldsymbol{u}_j \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{u}_j\right ) \,\textrm{d} V\\ &=-2 \sum _{j=1,2,3} \rho _j v_j \int _{V_j} \frac {\partial u_{i, k}}{\partial x_j} \frac {\partial u_{i, j}}{\partial x_k} \, \textrm{d} V \approx -2 \sum _{j=1,2,3} \rho _j v_j \int _{V_j}\left (\frac {\partial ^2 \varPhi _j}{\partial x_j \partial x_k}\right )^2 \, \textrm{d} V, \end{aligned} \end{equation}

where $\varPhi _j$ is the potential function of fluid $j$ . Thus, we obtain the internal damping

(3.18) \begin{align} \lambda ^{{Int}} & = \frac {\sum _{j=1}^3 \int {\rm \mu} _j\left (\frac {\partial ^2 \varPhi _j}{\partial x_j \partial x_k}\right )^2 \, \textrm{d} V_j}{\sum _{j=1}^3 \int \rho _j\left (\frac {\partial \varPhi _j}{\partial x_j}\right )^2 \, \textrm{d} V_j} \nonumber\\[4pt]& = \frac { {\begin{array}{c}2 k_{1,n}^2\Big (\epsilon ^2 \rho _1 v_1 \operatorname {coth}\left (h_1 k_{1,n}\right ) +\left (1+\epsilon ^2\right ) \rho _2 v_2 \operatorname {coth}\left (h_2 k_{1,n}\right ) \\ +\rho _3 v_3 \operatorname {coth}\left (h_3 k_{1,n}\right )-2 \epsilon \rho _2 v_2 \operatorname {csch}\left (h_2 k_{1,n}\right )\Big )\end{array}}}{ {\begin{array}{c}\epsilon ^2 \rho _1 \operatorname {coth}\left (h_1 k_{1,n}\right )+\left (1+\epsilon ^2\right ) \rho _2 \operatorname {coth}\left (h_2 k_{1,n}\right ) \\ +\rho _3 \operatorname {coth}\left (h_3 k_{1,n}\right )-2 \epsilon \rho _2 \operatorname {csch}\left (h_2 k_{1,n}\right )\end{array}}}, \nonumber\\\end{align}

with $\epsilon _{ \pm }=a_n \pm \sqrt {a_n^2+b_n}$ , $a_n=({1}/{2}) (1/\rho _{B,n}-(\rho _2-\rho _3)/(\rho _{A,n}(\rho _1-\rho _2)) )$ , $b_n=(\rho _2-\rho _3)/(\rho _1-\rho _2)$ . Since the frequency range investigated in this study closely matches the natural frequency corresponding to the (-) mode, we adopt the value of $\epsilon _-$ . Note that $\lambda ^{{BL}}$ and $\lambda ^{{IL}}$ have been determined by Herreman et al. (Reference Herreman, Wierzchalek, Horstmann, Cappanera and Nore2023). Here, $\lambda ^{{BL}} = \lambda _1^{{BL}}+\lambda _2^{{BL}}+\lambda _3^{{BL}}$ , where the contributions from the wall boundary layers are

(3.19) \begin{align} &\lambda _1^{{BL}} =\frac {\sqrt {2}}{4} \frac {1}{k_{1,n} R} \frac {\rho _1 \varOmega ^2 \sqrt {v_1 \varOmega }}{g\left (\rho _2-\rho _3+\epsilon ^2 \left (\rho _1-\rho _2\right )\right )} \nonumber\\[1.5pt] & \qquad\qquad\times\, \epsilon ^2\left [k_{1,n}\left (R-h_1\right ) \operatorname {csch}^2\left (k_{1,n} h_1\right )+\frac {k_{1,n}^2 R^2+1}{k_{1,n}^2 R^2-1} \operatorname {coth}\left (k_{1,n} h_1\right )\right ]\!, \nonumber\\[1.5pt] &\lambda _2^{{BL}} = -\frac {\sqrt {2}}{4} \frac {1}{k_{1, n} R} \frac {\rho _2 \varOmega ^2 \sqrt {v_2 \varOmega }}{g\left (\rho _2-\rho _3+\epsilon ^2\left (\rho _1-\rho _2\right )\right )} \nonumber\\[1.5pt] &\qquad\qquad \times\, \left (\epsilon ^2+1\right )\left (k_{1, n} h_2 \operatorname {csch}^2\left (k_{1, n} h_2\right )-\frac {k_{1, n}^2 R^2+1}{k_{1, n}^2 R^2-1} \operatorname {coth}\left (k_{1, n} h_2\right )\right ) \nonumber\\[1.5pt] &\qquad\qquad+\,\frac {\sqrt {2}}{2} \frac {1}{k_{1, n} R} \frac {\epsilon \rho _2 \varOmega ^2 \sqrt {v_2 \varOmega }}{g\left (\rho _2-\rho _3+\epsilon ^2\left (\rho _1-\rho _2\right )\right )} \nonumber\\[1.5pt] &\qquad\qquad \times\, \left (k_{1, n} h_2 \operatorname {coth}\left (k_{1, n} h_2\right ) \operatorname {csch}\left (k_{1, n} h_2\right )-\frac {k_{1, n}^2 R^2+1}{k_{1, n}^2 R^2-1} \operatorname {csch}\left (k_{1, n} h_2\right )\right )\!, \nonumber\\[1.5pt] &\lambda _3^{{BL}} = \frac {\sqrt {2}}{4} \frac {1}{k_{1,n} R} \frac {\rho _3 \varOmega ^2 \sqrt {v_3 \varOmega }}{g\left (\rho _2-\rho _3+\epsilon ^2\left (\rho _1-\rho _2\right )\right )} \nonumber\\[1.5pt] &\qquad\qquad \times\, \left [k_{1,n}\left (R-h_3\right ) \operatorname {csch}^2\left (k_{1,n} h_3\right )+\frac {k_{1,n}^2 R^2+1}{k_{1,n}^2 R^2-1} \operatorname {coth}\left (k_{1,n} h_3\right )\right ]\!. \end{align}

The sum damping of contributions from the interfacial boundary layers is

(3.20) \begin{equation} \lambda ^{{IL}} = \frac {\sqrt {2}}{4} \frac {\varOmega ^2 \sqrt {\varOmega }}{g(\rho _2-\rho _3+\epsilon ^2(\rho _1-\rho _2))} \sum _{j=1,2} \varLambda _{j j+1}^2\left(\frac {1}{\rho _j \sqrt {v_j}}+\frac {1}{\rho _{j+1} \sqrt {v_{j+1}}}\right)^{-1}, \end{equation}

with

(3.21) \begin{equation} \begin{gathered} \varLambda _{12} = \epsilon (\operatorname {coth}(k_{1,n} h_1)+\operatorname {coth}(k_{1,n} h_2))-\operatorname {csch}(k_{1,n} h_2),\\ \varLambda _{23}=(\operatorname {coth}(k_{1,n} h_2)+\operatorname {coth}(k_{1,n} h_3))-\epsilon \operatorname {csch}(k_{1,n} h_2). \end{gathered} \end{equation}

Usually, interface coupling effects exist in inertial terms. The strength of these coupling effects depends on $\rho _{A, n}$ and $\rho _{B, n}$ . However, it is necessary to make corrections for the interface coupling effects caused by viscosity, because the motion contributes to inertial and viscous forces. Therefore, the effect of interfacial coupling is also considered when viscosity correction terms with Herreman’s viscous decay rate $\lambda$ (Herreman et al. Reference Herreman, Wierzchalek, Horstmann, Cappanera and Nore2023) are added. The following modified equations for the three-layer fluid system have been derived:

(3.22a) \begin{align} & \underbrace {\frac {\textrm{d}^2 \beta _{A, n}}{\, \textrm{d} \phi ^2}-\rho _{A, n} \frac {\, \textrm{d}^2 \beta _{B, n}}{\, \textrm{d} \phi ^2}}_{\text{Inertial term}}+\underbrace {2 \varUpsilon \frac {\, \textrm{d} \beta _{A, n}}{\, \textrm{d} \phi }-2 \rho _{A, n} \varUpsilon \frac {\, \textrm{d} \beta _{B, n}}{\, \textrm{d} \phi }}_{\text{Viscous term}}+\underbrace {\frac {\varOmega _{A, n}^2}{\varOmega ^2} \beta _{A, n}}_{{\rm Restoring \,term}} \underbrace {-F_{A, n} {d}_s \sin \phi }_{\text{Excitation term}}=0, \end{align}

(3.22b) \begin{align} & \underbrace {\frac {\, \textrm{d}^2 \beta _{B, n}}{\, \textrm{d} \phi ^2}-\rho _{B, n} \frac {\, \textrm{d}^2 \beta _{A, n}}{\, \textrm{d} \phi ^2}}_{\text{Inertial term}}+\underbrace {2 \varUpsilon \frac {\, \textrm{d} \beta _{B, n}}{\, \textrm{d} \phi }-2 \rho _{B, n} \varUpsilon \frac {\, \textrm{d} \beta _{A, n}}{\, \textrm{d} \phi }}_{\text{Viscous term}}+\underbrace {\frac {\varOmega _{B, n}^2}{\varOmega ^2} \beta _{B, n}}_{{\rm Restoring \,term}} \underbrace {-F_{B, n} {d}_s \sin \phi }_{\text{Excitation term}}=0, \end{align}

where $\varUpsilon =\lambda /\varOmega$ describes the ratio of viscous force to inertial force.

Equation (3.22) are equilibrium equations that describe the dynamics of interfaces A and B. The two equations are coupled through interfacial coupling coefficients. The viscous contribution is introduced via a dissipation rate, representing energy loss due to viscosity. The restoring term accounts for the combined effects of gravity and interfacial tension, providing the ’elastic response’ of the interface, while the external term models the forcing induced by the orbital motion of the container. In two-layer systems with low viscosity, it is appropriate to model the viscous effect by adding a corrective term of the form $\varUpsilon \, \textrm{d} \beta _{A,n} / \textrm{d} \phi$ , ensuring that the energy dissipation rate matches $\varUpsilon$ . However, this single correction is insufficient in a three-layer system with finite viscosity. The interfacial coupling alters the decay characteristics, and the modified system no longer exhibits a decay rate equal to $\varUpsilon$ . To rectify this, an additional term $2 \rho _{A,n} \varUpsilon \, \textrm{d} \beta _{B,n} / \textrm{d} \phi$ is introduced, ensuring consistency with the prescribed damping rate. By solving (3.22), we can obtain an analytical expression for the wave amplitude of the oscillation in a three-layer fluid system according to potential flow theory (assuming counterclockwise). The corresponding expression for the interface fluctuation is given in (D1).

3.2. The semi-analytical model

The current modified model cannot satisfy the viscous boundary conditions at the interface and, therefore, cannot accurately capture the frequency shift. Bauer & Eidel (Reference Bauer and Eidel1999) and Ibrahim (Reference Ibrahim2005) established a semi-analytical model for single-layer fluid that strictly satisfies sidewall boundary conditions, but is inconvenient to extend to multi-layer fluid systems. Kumar & Tuckerman (Reference Kumar and Tuckerman1994) and Batson et al. (Reference Batson, Zoueshtiagh and Narayanan2013) applied the Floquet series to the Faraday wave problem analytically in an infinite domain. The approach neglected viscous shear at the sidewall while keeping the strict interface boundary conditions and no-slip conditions at the top and bottom walls. This simplification is adequate for most regions far from the sidewall boundary layer. Experiments in cylindrical containers confirmed the validity of this theoretical approach (Batson et al. Reference Batson, Zoueshtiagh and Narayanan2013; Ward et al. Reference Ward, Zoueshtiagh and Narayanan2019c ).

Considering the $O(\varepsilon )$ -order approximation, $\xi _A \approx h_1+\eta _A$ , $\xi _B \approx h_1+h_2+\eta _B$ , the linear solutions can be expressed in the form of a Fourier series expansion

(3.23) \begin{align} \bar {w}_j(r, \phi , z) &= \sum _{n=1}^{\infty } W_{j, n}(z) J_1\left (k_{1, n} r\right ) e^{i \phi }, \,\, \bar {p}_j(r, \phi , z)=\sum _{n=1}^{\infty } P_{j, n}(z) J_1\left (k_{1, n} r\right ) e^{i \phi }, \nonumber\\ \bar {\eta }_A(r, \phi ) &= \sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right ) e^{i \phi }, \,\, \bar {\eta }_B(r, \phi )=\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right ) e^{i \phi }, \end{align}

with $\phi =\varOmega t-\varphi$ . In addition, $\bar {w}_j$ is the velocity in the $z$ -direction and $\bar {p}_j$ is the pressure of fluid $j$ in the form of complex functions. In addition, $\varepsilon _{1,n}$ represents the nth zero of the first-order Bessel function derivative ${J}^{\prime}_1 ( \varepsilon _{1,n} )$ , $k_{1,n}=\varepsilon _{1,n}/R$ , because the slip boundary on the side wall has been satisfied. Applying the divergence operator to the Navier–Stokes equations (3.2b ) and simplifying, we obtain

(3.24) \begin{equation} \frac {\textrm{d}^2 P_{j, n}}{\textrm{d} z^2}-k^2_{1,n}P_{j, n}=0. \end{equation}

The general solution of the equation can be solved,

(3.25) \begin{equation} P_{j, n}=C_{j, n} e^{-k_{1, n} z}+D_{j, n} e^{k_{1, n} z^2}, \end{equation}

where $C_{j, n}$ and $D_{j, n}$ are undetermined coefficients.

Substituting the series solution (3.23) into (3.2b ), we obtain

(3.26) \begin{equation} \frac {\textrm{d}^2 W_{j, n}}{\textrm{d} z^2}-\left (k_{1, n}^2+\frac {\textrm{i} \varOmega }{v_j}\right ) W_{j, n}=\frac {1}{{\rm \mu} _j} \frac {\textrm{d} P_{j, n}}{\textrm{d} z}. \end{equation}

Solving (3.26), we obtain the general solution for $W_{j, n}$ ,

(3.27) \begin{equation} W_{j, n}=A_{j, n} \textrm{e}^{-\sqrt {k_{1, n}^2+\textrm{i} \frac {\varOmega }{v_j} }z}+B_{j, n} \textrm{e}^{\sqrt {k_{1, n}^2+\textrm{i} \frac {\varOmega }{v_j} }z}+C_{j, n} \frac {-\textrm{i} k_{1, n}}{\rho _j \varOmega } \textrm{e}^{-k_{1, n} z}+D_{j, n} \frac {\textrm{i} k_{1, n}}{\rho _j \varOmega } \textrm{e}^{k_{1, n} z}, \end{equation}

where $A_{j, n}$ and $B_{j, n}$ are undetermined coefficients.

The boundary conditions are also expressed as follows.

1. (i) Velocity boundary conditions:

   (3.28) \begin{align} \left .W_{1, n}\right |_{z=0} &= 0,\left .\frac {\textrm{d} W_{1, n}}{\textrm{d} z}\right |_{z=0}=0, \nonumber\\ \left .W_{1, n}\right |_{z=h_1}&= \left .W_{2, n}\right |_{z=h_1}=\textrm{i} \varOmega \zeta _{A, n},\left .\frac {\textrm{d} W_{1, n}}{\textrm{d} z}\right |_{z=h_1}=\left .\frac {\textrm{d} W_{2, n}}{\textrm{d} z}\right |_{z=h_1}, \nonumber\\[4pt] \left .W_{1, n}\right |_{z=h_1+h_2} &= \left .W_{2, n}\right |_{z=h_1+h_2}=\textrm{i} \varOmega \zeta _{B, l},\left .\frac {\textrm{d} W_{1, n}}{\textrm{d} z}\right |_{z=h_1+h_2}=\left .\frac {\textrm{d} W_{2, n}}{\textrm{d} z}\right |_{z=h_1+h_2}, \nonumber\\ \left .W_{3, n}\right |_{z=\sum h_i} &= 0,\left .\frac {\textrm{d} W_{3, n}}{\textrm{d} z}\right |_{z=\sum h_i}=0 . \end{align}

2. (ii) Shear stress boundary conditions:

   (3.29) \begin{align} \left .{\rm \mu} _1\left (k_{1, n}^2 W_{1, n}+\frac {\textrm{d}^2 W_{1, n}}{\textrm{d} z^2}\right ) \right |_{z=h_1}-\left .{\rm \mu} _2\left (k_{1, n}^2 W_{2, n}+\frac {\textrm{d}^2 W_{2, n}}{\textrm{d} z^2}\right ) \right |_{z=h_1} &= 0, \nonumber\\[6pt] \left .{\rm \mu} _2\left (k_{1, n}^2 W_{2, n}+\frac {\textrm{d}^2 W_{2, n}}{\textrm{d} z^2}\right ) \right |_{z=h_1+h_2}-\left .{\rm \mu} _3\left (k_{1, n}^2 W_{3, n}+\frac {\textrm{d}^2 W_{3, n}}{\textrm{d} z^2}\right ) \right |_{z=h_1+h_2} &= 0 . \end{align}

3. (iii) Normal stress boundary conditions:

   (3.30) \begin{equation} \begin{aligned} & \left .P_{1, n}\right |_{z=h_1}-\left .P_{2, n}\right |_{z=h_1}-\left .2\left ({\rm \mu} _1-{\rm \mu} _2\right ) \frac {\textrm{d} W_{1, n}}{\textrm{d} z}\right |_{z=h_1} \\[3pt]& -\big [(\rho _1-\rho _2) g+\sigma _A k_{1, n}^2\big ] \zeta _{A, n}=-\frac {R\left (\rho _1-\rho _2\right ) {d}_s \varOmega ^2}{\big (k_{1, n}^2 R^2-1\big ) J_1\big (k_{1, n} R\big )} \\[3pt] & \left .P_{2, n}\right |_{z=h_1+h_2}-\left .P_{3, n}\right |_{z=h_1+h_2}-\left .2\left ({\rm \mu} _2-{\rm \mu} _3\right ) \frac {\textrm{d} W_{3, n}}{\textrm{d} z}\right |_{z=h_1+h_2} \\[3pt] & -\big [(\rho _2-\rho _3) g+\sigma _B k_{1, n}^2\big ] \zeta _{B, l}=-\frac {R\left (\rho _2-\rho _3\right ) {d}_s \varOmega ^2}{\big (k_{1, n}^2 R^2-1\big ) J_1\big (k_{1, n} R\big )}. \end{aligned} \end{equation}

Substituting (3.25), (3.27) into (3.28), (3.29), (3.30) allows us to solve for $A_{j, n}$ , $B_{j, n}$ , $C_{j, n}$ , $D_{j, n}$ , $\zeta _{A, n}$ and $\zeta _{B, n}$ , where $j=1,2,3$ . The boundary conditions (3.28), (3.29), (3.30) provide 14 equations, allowing the determination of 14 unknowns for any given positive integer $n$ . The counterclockwise solutions of $\eta _A$ and $\eta _B$ are the real solution,

(3.31) \begin{eqnarray} \eta _A &=&\delta _A \cos \left (\phi +\varphi _A\right )=\textrm{Re}\left (\bar {\eta }_A\right ) \nonumber\\ &=& \left (\textrm{Re}\left (\sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right )\right ) \cos \phi -\textrm{Im}\left (\sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right )\right ) \sin \phi \right ), \nonumber\\ \eta _B &=& \delta _B \cos \left (\phi +\varphi _B\right )=\textrm{Re}\left (\bar {\eta }_B\right )\nonumber\\ &=& \left (\textrm{Re}\left (\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right )\right ) \cos \phi -\textrm{Im}\left (\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right )\right ) \sin \phi \right ), \end{eqnarray}

with

(3.32) \begin{equation} \begin{gathered} \delta _A(r)=\left |\sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right )\right |, \delta _B(r)=\left |\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right )\right |, \\[4pt] \varphi _A=\arctan \frac {\textrm{Im}\left (\sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right )\right )}{\textrm{Re}\left (\sum _{n=1}^{\infty } \zeta _{A, n} J_1\left (k_{1, n} r\right )\right )}, \varphi _B=\arctan \frac {\textrm{Im}\left (\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right )\right )}{\textrm{Re}\left (\sum _{n=1}^{\infty } \zeta _{B, n} J_1\left (k_{1, n} r\right )\right )}. \end{gathered} \end{equation}

4. Scaling laws

Figure 3.

Variation of amplitude ratio and phase difference with $Fr$ for different Reynolds numbers. Here, ${Re}_{j,{Air\mid P\mid S}}$ represents the Reynolds number of the $j$ th fluid layer in the Air∣P∣S fluid system. At $Fr = 0.01$ , the Reynolds numbers for the three fluid layers are ${Re}_1 = 90$ , ${Re}_2 = 112$ and ${Re}_3 = 303$ . When $Fr=0.1$ , they increase to ${Re}_1 = 283$ , ${Re}_2 = 354$ and ${Re}_3 = 957$ , respectively, indicating that ${Re}_i \gg 1$ in all cases. Other parameters are the same as those of Air∣P∣S.

We define the Froude number ${\textit{Fr}} = {d}_s \varOmega ^2 / (2g)$ , which represents the ratio of inertial force to gravitational force, and Reynold number ${\textit{Re}}_j = \varOmega R^2 / \nu _{\!j}$ , which represents the ratio of inertial force to viscous force. With the Air∣P∣S system as a reference, we calculate the variation of the amplitude ratio $\delta _A / \delta _B$ and phase difference $\varphi / 2\pi$ with $Fr$ for different ${\textit{Re}}_j$ . As shown in figure 3, the amplitude ratio and phase difference reach their maximum values at resonance from the calculations of $\S$ 3.2. At resonance, the amplitude ratio $\delta _A / \delta _B$ and the phase difference $\varphi / 2\pi$ increase with increasing ${\textit{Re}}_j$ . When ${\textit{Re}}_j = \infty$ , the wave amplitude ratio exhibits non-physical singularities and the phase difference undergoes an abrupt shift at resonance. However, the viscous model eliminates these singularities in amplitude and discontinuities in phase difference, providing a more realistic description of near resonance.

While the semi-analytical model offers higher accuracy by fully resolving viscous effects, its analytical complexity limits its utility for scaling analysis. The modified model, however, enables a convenient perturbative approach, making it well suited for extracting the scaling behaviours of the system. We separately decompose $\beta _{A,n}$ and $\beta _{B,n}$ into components $\beta _{A,n} =\beta _{A,S,n}\sin \phi +\beta _{A,C,n}\cos \phi$ and $\beta _{B,n} =\beta _{B,S,n}\sin \phi +\beta _{B,C,n}\cos \phi$ . However, we consider the characteristics of solutions far from resonance. The small perturbation expansion of (D1) is conducted with $\varUpsilon =\lambda /\varOmega \ll 1$ . From (D5), which is a non-resonant asymptotic solution, we know $\beta _{A,C,n},\beta _{B,C,n} \sim \varUpsilon$ . Therefore, $\eta _A,\eta _B \sim \varUpsilon$ far from resonance. However, we apply the small perturbation expansion at resonance satisfying $\gamma _n=0$ . From (D6), which is the non-resonant asymptotic solution, we know $\beta _{A,C,n},\beta _{B,C,n} \sim \varUpsilon ^{-1}$ . Therefore, $\eta _A,\eta _B \sim \varUpsilon ^{-1}$ at resonance.

At resonance, when $\left | \tan \varphi _A \right | \ll 1$ and $\left |\tan \varphi _B \right | \ll 1$ , the phase difference reaches an extremum $\varphi _{\textit{max}}$ , which approximately satisfies

(4.1) \begin{equation} \pi - \Delta \varphi _{\textit{max}} \approx \frac {\beta _{A, S, n}}{\beta _{A, C, n}}-\frac {\beta _{B, S, n}}{\beta _{B, C, n}}, \end{equation}

because $\Delta \varphi _{\textit{max}}$ is close to $\pi$ , such that $\Delta \varphi _{\textit{max}} \approx \tan \Delta \varphi _{\textit{max}} \approx \tan \varphi _B - \tan \varphi _A$ . Moreover, since the resonant wave amplitudes scale as $\eta _A, \eta _B \sim \varUpsilon ^{-1}$ , the resulting phase difference follows the scaling law $\pi -\Delta \varphi _{\textit{max}} \sim \varUpsilon$ .

Using the semi-analytical model in $\S$ 3.2, we compute $\delta _A$ and $\Delta \phi$ over the range $Fr = 0\,{\sim}\,0.1$ , from which the first-order resonance peak values $\delta _{A,{\textit{max}}}$ and $\Delta \phi _{{\textit{max}}}$ are extracted. Figure 4 displays the relationships between the first-order resonant amplitude $\delta _{A,\textit{max}}$ of interface $A$ , the maximum phase difference $\Delta \varphi _{{\textit{max}}}$ , and the viscous damping $\varUpsilon$ of the fluid system. In figure 4(a), the first-order resonant amplitude of interface $A$ increases with the rise of $\varUpsilon ^{-1}$ , indicating that higher viscous damping in the fluid system leads to smaller resonant amplitude. When the fluid system exhibits low viscous damping $\varUpsilon ^{-1}\gg 1$ , the first-order resonant amplitude satisfies $\delta _{A,{\textit{max}}}/{\textrm d}s\sim \varUpsilon ^{-1}$ , as the condition $\eta_A, \eta_B \sim \varUpsilon ^{-1}$ is satisfied at resonance. This conclusion is consistent with Horstmann et al. (Reference Horstmann, Anders, Kelley and Weier2021), namely $\delta _{\textit{max}}/R \sim \sqrt {{\textit{Re}}_j}$ , because $\lambda \sim {\textit{Re}}_j^{-1/2}$ . The maximum phase difference between interface $A$ and interface $B$ approaches $\pi$ as $\varUpsilon$ decreases. Under low viscous damping conditions ( $\varUpsilon \ll 1$ ), the maximum phase difference satisfies $\pi -\Delta \varphi _{{\textit{max}}}\sim \varUpsilon$ , as illustrated in figure 4(b). Additionally, as the coupling coefficients $\rho _{A,n}$ and $\rho _{B,n}$ increase, the maximum amplitude of interface $A$ decreases, and the maximum phase difference also decreases. This behaviour is related to energy transfer between interfaces A and B.

Figure 4.

Normalised resonance amplitude and phase difference maximum. (a) Dimensionless resonance amplitude at interface $A$ as a function of $\varUpsilon ^{-1}$ . The red arrows indicate the direction of enhanced coupling. (b) Maximum of phase difference $\Delta \varphi _{{\textit{max}}}$ as a function of $\varUpsilon$ . The viscosity ${\textit{Re}}_j=\varOmega R^2/\nu _{\!j}$ ranges from 0.01 to 10 times the original fluid viscosity of Air∣P∣S. Here, ${\textit{Re}}_{j,{Air\mid P\mid S}}$ represents the Reynolds number of the $j$ th fluid layer in the Air∣P∣S fluid system. Other parameters are the same as those of Air∣P∣S.

Figure 5.

Impact of additional fluid layers on interfacial kinetic energy. (a) Comparison of interfacial kinetic energy as a function of frequency in a three-layer fluid system (Air∣P∣S system) and a two-layer fluid system (Air∣P system) with 50 cst of silicone oil and $h_2$ = 20 mm. (b) Relationship between the kinetic energy growth rate and frequency for different Reynolds numbers ${\textit{Re}}_j$ . (c) Coupling parameters $\rho _A,\rho _B$ as a function of $h_2/h_{2,{Air\mid P\mid S}}$ . (d) Relationship between the maximum kinetic energy growth rate and $\varUpsilon ^{-1}$ . ${\textit{Re}}_{j,{Air\mid P\mid S}}$ represents the Reynolds number of the $j$ th fluid layer in the Air∣P∣S fluid system. Other parameters are the same as those of Air∣P∣S.

Therefore, we further explore the changes in kinetic energy caused by interface coupling, using $\boldsymbol{\tilde {u}}^2_B=\boldsymbol{u}^2_B/({d}_s\varOmega )^2$ to quantify the interfacial kinetic energy at interface $B$ , where $\boldsymbol{\tilde {u}}_B$ represents the interfacial normalised velocity. Furthermore, dimensional analysis suggests that the interfacial velocity $\boldsymbol{u_B}^2$ satisfies the scaling relation $|\boldsymbol{u}_B| \sim \varOmega \eta _B$ . Therefore, $\boldsymbol{u_B}^2$ satisfies $\boldsymbol{u}_B^2 \sim \eta _B^2 \sim \varUpsilon ^{-2}$ following the resonance condition. Compared with the two-layer fluid system (Air∣P), introducing an additional fluid interface in the three-layer configuration (Air∣P∣S) significantly changes the interfacial dynamics. This additional P–S interface is a conduit for enhanced energy exchange between layers. As illustrated in figure 5(a), the presence of the intermediate silicone oil (S) layer leads to a pronounced jump in kinetic energy at interface $B$ (Air–Paraffin oil interface) when the system is driven at the resonant frequency $\varOmega = \varOmega _A$ , corresponding to the natural frequency of the Air∣P configuration. This abrupt increase is reminiscent of the classic dead water phenomenon first described by Ekman (Reference Ekman1904) in his doctoral thesis, wherein a vessel moving through a stratified fluid experiences a significantly greater drag than in an homogeneous fluid (Mercier, Vasseur & Dauxois Reference Mercier, Vasseur and Dauxois2011). The enhanced resistance arises from wave generation at the fluid interface, which extracts energy from the moving body. The effect is most pronounced when the vessel speed is below the maximum phase speed of interfacial waves, with the interfacial kinetic energy $\boldsymbol{u}^2$ as a proxy for the drag (Fourdrinoy et al. Reference Fourdrinoy, Dambrine, Petcu, Pierre and Rousseaux2020). Figure 5(b) reveals that this kinetic energy jump becomes more abrupt as viscosity decreases, indicating that viscous shear plays a mitigating role. In the low-viscosity limit ( $\varUpsilon ^{-1} \gg 1$ ), the energy growth rate follows $\tilde {\boldsymbol{u}}^2_{B,3-\textit{layer,max}} / \tilde {\boldsymbol{u}}^2_{B,2-\textit{layer}} \sim \varUpsilon ^{-2}$ , consistent with $\boldsymbol{u}^2_B \sim \varUpsilon ^{-2}$ . Moreover, as shown in figure 5(c, d), stronger interfacial coupling amplifies this resonant energy growth, particularly when the intermediate fluid layer is thin. Introducing a third, thin interfacial layer substantially increases energy transfer and drag by exciting interfacial waves at multiple boundaries. This enhanced coupling effectively mimics the sharp density transition in classical two-layer systems, reinforcing the analogy to Ekman’s dead water.

Figure 6.

BOS images corresponding to the different viscosity of fluid 1 and shaking frequency, where $\nu _{1,\textit{P}\mid \textit{S}\mid \textit{G}}$ is the kinematic viscosity of glycerol. Other parameters are the same as those of P∣S∣G.

Through calculations using the theoretical model, we observed BOS patterns with spiral structures similar to those found in experiments on two-layer fluid systems (Horstmann et al. Reference Horstmann, Anders, Kelley and Weier2021), as shown in figure 6. The number of spiral arms increases with the sloshing frequency. Additionally, as the viscosity of the bottom fluid layer increases, the system gradually transitions from underdamped to overdamped characteristics.

5. Comparison with experimental data

To measure the amplitude and phase quantitatively of the interfaces, we use ultrasonic techniques previously introduced by Horstmann et al. (Reference Horstmann, Wylega and Weier2019) to measure the wave amplitudes of a single fluid interface. This approach enables the simultaneous measurement of fluctuations in the two fluid interfaces. For ultrasonic measurement technology, we use the Air∣P∣S fluid system to avoid extreme damping of ultrasonic waves in the fluid. In figures 7(a)–7(c) and figures 8(a)–8(c), the kinematic viscosity of silicone oil is 50 cst. In figures 7(d)–7(f) and figures 8(d)–8(f), the kinematic viscosity of silicone oil is 500 cst. Additionally, the initial thickness of the paraffin oil layer is 20 mm in figures 7(a, d) and 8(a, d), 25 mm in figures 7(b, e) and 8(b, e), and 30 mm in figures 7(c, f) and 8(c, f). The ultrasonic probe is fixed at $r_0 = 30\,{\textrm{mm}}$ .

Figure 7.

Comparison of experimental measurements and theoretical predictions of the wave amplitudes, $\delta _A$ and $\delta _B$ , in the three-layer fluid system. We select three typical thicknesses of middle-layer paraffin oil, (a, d) 20 mm, (b, e) 25 mm and (c ,f) 30 mm. The viscosity of silicone oil is 50 cst in figures 7(a)–7(c) and 500 cst in figures 7(d)–7(f). The hollow circles denote experimental data, the dashed lines show the data predicted by the modified model and the solid lines represent the predictions calculated by the semi-analytic model.

Figure 8.

Comparison of experimental measurements and theoretical predictions of the phase difference $\Delta \varphi$ of interface waves in the three-layer fluid system. We select three typical thicknesses of middle-layer paraffin oil, (a, d) 20 mm, (b, e) 25 mm and (c, f) 30 mm. The viscosity of silicone oil is 50 cst in figures 8(a)–8(c) and 500 cst in figures 8(d)–8(f). $\varphi _{A}$ and $\varphi _{B}$ represent the phase differences between interfaces and the shaker. The red ( $\varphi _{A}$ ) and blue ( $\varphi _{B}$ ) dashed lines in the graph represent the phase differences between interfaces and the shaker, as calculated by the semi-analytical model.

Figure 9.

Comparison of experimental measurements and theoretical predictions of the resonance frequency of interface $A$ .

We perform a comparative evaluation of the model predictions against experimental data for the three-layer fluid system, as plotted in figures 7 and 8. The amplitudes at $r=r_0$ on both interfaces are measured using ultrasonic techniques, as indicated in figure 7 by the hollow circles. The results calculated using the semi-analytical and modified models agree well with the experimental data. Similar to the two-layer fluid system, the semi-analytical model captures more accurately the frequency shift caused by viscosity. It confirms that the primary boundary conditions are associated with the interface in the flow field near the interface and away from the sidewall boundary layer. The influence of the convection and other higher-order terms drives the deviation from experimental measurements and the two models. On the one hand, the premise of the model in this article is based on the linear wave assumption; therefore, when the amplitude is large, there may be some deviation from the experiments. On the other hand, the amplitude of resonance may lead to a smaller amplitude than linear waves due to superharmonic resonance (Bongarzone et al. Reference Bongarzone, Guido and Gallaire2022). Due to viscosity, phase differences, $\varphi _{A}$ and $\varphi _{B}$ , exist between the two interfaces and the shaker, as illustrated in figure 1(b). As shown in figure 8, $\Delta \varphi =\varphi _{B}-\varphi _{A}$ is well captured by ultrasound techniques. Both the semi-analytical model and the modified model can capture the phase difference. However, the results of the semi-analytical model are closer to the experimental values, attributed to fully satisfying the fluid dynamics equations and interface boundary conditions. Figure 8 shows that interface $A$ primarily contributes to the phase difference between the two interfaces. Passing the first resonance point, interface $A$ introduces a phase shift close to $\pi /2$ (exactly $\pi /2$ in the inviscid case) relative to the sloshing platform, while interface $B$ remains far from its resonance point. To highlight the advantage of the semi-analytical model in capturing the viscous frequency shift, figure 9 compares its predictions, those of the modified model and experimental measurements. The semi-analytical results show significantly better agreement with the experimental data.

Figure 10.

Normalisation of the sum of gradients of internal waves at two interfaces. (a) Experimental results were obtained by using the BOS method with $\varOmega = 100\,\textrm{rpm}$ . (b) Reconstructed results from the semi-analytical model with $\varOmega = 100\,\textrm{rpm}$ . (c) Reconstructed results from the modified model with $\varOmega = 100\,\textrm{rpm}$ .

By numerically solving the semi-analytical model and the modified model, the spiral arm patterns are observed in the three-layer fluid system, as shown in figure 6. To validate this theoretical finding, we applied the BOS technique to capture wave gradient information in a two-dimensional plane for the P∣S∣G fluid system. The results in figure 10 show that both the semi-analytical and modified models provide reliable predictions not only for the wave height, but also for the estimation of the interface except near the sidewalls from figure 10. We experimentally measure the gradient distribution of the interface except near the sidewalls. Using the BOS technique, we experimentally measured the gradient distribution of interface waves, as shown in figure 10(a). Our measured data validate the spiral arm structure predicted by both the semi-analytical model (figure 10 b) and the modified model (figure 10 c). Moreover, the results of the semi-analytic model are more consistent with the measured results, as evident in the angle of the spiral arm caused by the viscosity. The discrepancy can be attributed to the effects of BOS and the normalisation process. Specifically, local defects in the container or the speckle paper may result in excessive local measurement gradients due to the BOS, which impacts the alignment between the normalised values and the theoretical predictions. However, it is important to note that this issue does not compromise the overall validity of the comparative analysis.

Both theoretical models presented in this study rely on the assumption of linear waves, i.e. the wave amplitude is much smaller than the wave height and wavelength. The modified model, which introduces a viscous damping term into the inviscid solution, is valid under the condition ${\textit{Re}}_j \gg 1$ . The non-dimensional parameter $\varUpsilon$ can quantify this requirement more precisely, with $\varUpsilon \ll 1$ equivalent to ${\textit{Re}}_j \gg 1$ . Figure 4 of the manuscript shows that the modified model remains accurate when $\varUpsilon \lt 0.01$ . In contrast, the semi-analytical model does not rely on the assumption ${\textit{Re}}_j \gg 1$ , as it directly solves the linearised Navier–Stokes equations with viscous terms retained. However, this model is constrained by the boundary conditions at the lateral wall. Since the solution is constructed using a Fourier series expansion in an azimuthally unbounded domain, it effectively captures the dynamics in the central region far from the container sidewall. This is also why ultrasonic probes measure the interfacial amplitude away from the boundary layer instead of tracking the contact line. Additional boundary conditions must be imposed on both models if a pinned contact line is present. In such cases, it is often necessary to solve an integro-differential system, as shown by Shao et al. (Reference Shao, Wilson, Saylor and Bostwick2021). Moreover, both models are first-order approximations, neglecting mean drift velocities arising from nonlinear convective terms (Bongarzone & Gallaire Reference Bongarzone and Gallaire2024). Capturing such effects would require extending the expansion to the second order. At least third-order approximations may be necessary for pronounced superharmonic responses with multiple wave crests (Bongarzone et al. Reference Bongarzone, Guido and Gallaire2022).