2.1 Statement of the problem

Figure 1 shows the set-up for our theoretical framework, analogous to the description given by Reclari (Reference Reclari2013). We define an ideal cylinder of radius $R$, containing two liquid layers (subscripts $i=1,2$) of heights $h_{1},h_{2}$, kinematic viscosities $\unicode[STIX]{x1D708}_{1},\unicode[STIX]{x1D708}_{2}$ and densities $\unicode[STIX]{x1D70C}_{1},\unicode[STIX]{x1D70C}_{2}$, where $\unicode[STIX]{x1D70C}_{1}<\unicode[STIX]{x1D70C}_{2}$ must be fulfilled to ensure a stable vertical stratification. Gravity $\boldsymbol{g}$ acts in the negative $z$-direction. In cylindrical coordinates ($r,\unicode[STIX]{x1D711},z$) the interface is located at $z=\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D711},t)$. Interfacial tension $\unicode[STIX]{x1D6FE}$ is considered along the liquid–liquid interface but capillary forces acting on the contact line are neglected; we assume that the interface may slide freely along the cylinder wall while maintaining a static contact angle of $90^{\circ }$ (no meniscus). The coordinate origin of the moving cylinder, the moving frame of reference $O(r,\unicode[STIX]{x1D711},z)$, is defined in the centre of the unperturbed interface $\unicode[STIX]{x1D702}$. In addition, we define the inertial frame of reference $O^{\prime }(r^{\prime },\unicode[STIX]{x1D711}^{\prime },z^{\prime })$ to describe the circulatory trajectories of the shaking table. The shaker is horizontally translated with a constant angular frequency $\unicode[STIX]{x1D6FA}$ along a circular path $\unicode[STIX]{x1D711}^{\prime }(t)=\unicode[STIX]{x1D6FA}t$ of shaking diameter $d_{s}$.

This formulation is characterised by eleven physical variables and three physical dimensions. Following the Buckingham $\unicode[STIX]{x03C0}$ theorem, the system can be uniquely described by eight independent dimensionless parameters. We define the following set of dimensionless numbers for our analysis:

(2.1)$$\begin{eqnarray}Fr={\displaystyle \frac{d_{s}\unicode[STIX]{x1D6FA}^{2}}{2g}},\quad E={\displaystyle \frac{d_{s}}{2R}},\quad Re_{i}={\displaystyle \frac{\unicode[STIX]{x1D6FA}R^{2}}{\unicode[STIX]{x1D708}_{i}}},\quad H_{i}={\displaystyle \frac{h_{i}}{R}},\quad Bo={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gR^{2}}{\unicode[STIX]{x1D6FE}}},\quad A={\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2}}}.\end{eqnarray}$$

The Froude number $Fr$ describes the forcing expressed by the ratio of inertial force and gravitational acceleration. The normalised shaking diameter determines the eccentricity $E$, while $Re_{i}$ are the phase-dependent Reynolds numbers, here weighting the cylinder radius with the boundary layer thicknesses $\unicode[STIX]{x1D6FF}_{i}=\sqrt{\unicode[STIX]{x1D708}_{i}/\unicode[STIX]{x1D6FA}}$. The importance of gravitational forces compared with interfacial tension forces is specified by the Bond number $Bo$. Finally, $H_{i}$ and $A$ are the geometrical aspect ratios and the Atwood number, respectively. It is the aim of our following analysis to relate these numbers to the resulting wave amplitudes and phase shifts.

2.2 Modal equations and solutions

The velocity vector $\boldsymbol{V}_{0}(t)$ describing the orbital container motion can be expressed in terms of cylindrical unit base vectors in the inertial frame of reference $O^{\prime }(r^{\prime },\unicode[STIX]{x1D711}^{\prime },z^{\prime })$

(2.2)$$\begin{eqnarray}\boldsymbol{V}_{0}(t)=\left(\begin{array}{@{}c@{}}-{\displaystyle \frac{d_{s}\unicode[STIX]{x1D6FA}}{2}}\sin (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711}^{\prime })\boldsymbol{e}_{r^{\prime }}\\[10.0pt] {\displaystyle \frac{d_{s}\unicode[STIX]{x1D6FA}}{2}}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711}^{\prime })\boldsymbol{e}_{\unicode[STIX]{x1D711}^{\prime }}\\ \end{array}\right).\end{eqnarray}$$

Due to this orbital motion the liquids in the cylinder constantly experience rotary acceleration in the moving frame of reference $O$ resulting in centripetal forces, which can be formulated as external volume forces in the equations of motion for the container-fixed coordinate system. The incompressible Euler equations for $O$ under the influence of gravity $\boldsymbol{g}$ and subject to $\boldsymbol{V}_{0}$ are

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{i}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{1}{2}}\unicode[STIX]{x1D735}|\boldsymbol{u}_{i}^{2}|-\boldsymbol{u}_{i}\times (\unicode[STIX]{x1D735}\times \boldsymbol{u}_{i})+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{i}}}\unicode[STIX]{x1D735}P_{i}-\boldsymbol{g}+\dot{\boldsymbol{V}}_{0}=0, & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{i}=0, & \displaystyle\end{eqnarray}$$

(compare with Kochin, Kibel & Roze (Reference Kochin, Kibel and Roze1964, chap. 2)), where $\boldsymbol{u}_{i}$ are the relative velocity fields and $P_{i}$ the pressures associated with both fluids $(i=1,2)$. By postulating irrotationality ($\unicode[STIX]{x1D735}\times \boldsymbol{u}_{i}=0$), the potential flow approximation can be applied stating that the flow fields $\boldsymbol{u}_{i}$ can be uniquely expressed as gradients of scalar potentials $\unicode[STIX]{x1D719}_{i}$

(2.5)$$\begin{eqnarray}\boldsymbol{u}_{i}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i}.\end{eqnarray}$$

By inserting (2.5) into (2.3) and (2.4) and integrating over the cylinder volume, a complete set of governing equations and boundary conditions for the scalar flow potentials $\unicode[STIX]{x1D719}_{i}$, the pressures $P_{i}$ and the interface position $\unicode[STIX]{x1D702}$ can be formulated

(2.6a)$$\begin{eqnarray}\displaystyle \text{Flow fields:}\!\quad & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}t}}-{\displaystyle \frac{d_{s}\unicode[STIX]{x1D6FA}^{2}r}{2}}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})+{\displaystyle \frac{1}{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i})^{2}+{\displaystyle \frac{P_{i}}{\unicode[STIX]{x1D70C}_{i}}}+gz=c_{i}(t),\end{eqnarray}$$

(2.6b)$$\begin{eqnarray}\displaystyle \displaystyle \text{Flow fields:}\!\quad & & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{i}=0,\end{eqnarray}$$

(2.6c)$$\begin{eqnarray}\displaystyle \displaystyle \text{Top wall:}\!\quad & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}=0|_{z=h_{1}},\end{eqnarray}$$

(2.6d)$$\begin{eqnarray}\displaystyle \displaystyle \text{Bottom wall:}\!\quad & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}=0|_{z=-h_{2}},\end{eqnarray}$$

(2.6e)$$\begin{eqnarray}\displaystyle \displaystyle \text{Sidewall:}\!\quad & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}r}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}r}}=0|_{r=R},\end{eqnarray}$$

(2.6f)$$\begin{eqnarray}\displaystyle \displaystyle \text{Interface:}\!\quad & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|_{z=\unicode[STIX]{x1D702}},\end{eqnarray}$$

(2.6g)$$\begin{eqnarray}\displaystyle \displaystyle \text{Interface:}\!\quad & & \displaystyle \unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}\boldsymbol{\cdot }{\displaystyle \frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}}{\sqrt{1+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}|^{2}}}}=P_{1}-P_{2}|_{z=\unicode[STIX]{x1D702}}.\end{eqnarray}$$

This formulation represents a class of spectral sloshing problems well known in the literature. Faltinsen & Timokha (Reference Faltinsen and Timokha2009) introduce different modal theories to find approximate analytical solutions for free-surface sloshing problems, mathematically quite similar to our two-layer formulation. In this study, we want to focus on first-order linear solutions, which are reasonable approximations as long as the wave amplitudes remain sufficiently small. In the linear approximation, the flow potentials $\unicode[STIX]{x1D719}_{i}$ can be described as superpositions of an infinite number of independent sloshing modes $m\in \mathbb{N}_{0}$ and $n\in \mathbb{N}_{1}$

(2.7a)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D719}_{1}(r,\unicode[STIX]{x1D711},z,t)=-\mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\unicode[STIX]{x1D6F7}_{mn}(\unicode[STIX]{x1D711},t){\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(z-h_{1})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)}}J_{m}\left(\unicode[STIX]{x1D716}_{mn}\frac{r}{R}\right),\end{eqnarray}$$

(2.7b)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D719}_{2}(r,\unicode[STIX]{x1D711},z,t)=\mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\unicode[STIX]{x1D6F7}_{mn}(\unicode[STIX]{x1D711},t){\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(z+h_{2})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}J_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right),\end{eqnarray}$$

(2.7c)$$\begin{eqnarray}\displaystyle & & \displaystyle \quad \text{with }\unicode[STIX]{x1D6F7}_{mn}(\unicode[STIX]{x1D711},t)=\unicode[STIX]{x1D6FC}_{mn}(t)\cos (m\unicode[STIX]{x1D711})+\unicode[STIX]{x1D6FD}_{mn}(t)\sin (m\unicode[STIX]{x1D711}),\end{eqnarray}$$

$\unicode[STIX]{x1D6FC}_{mn}(t)$

$\unicode[STIX]{x1D6FD}_{mn}(t)$

$J_{m}$

$m$

$\unicode[STIX]{x1D716}_{mn}$

$n$

$m$

$J_{m}^{^{\prime }}(\unicode[STIX]{x1D716}_{mn})=0$

$\unicode[STIX]{x1D716}_{mn}$

$m$

In order for the flow potentials (2.7) to fulfil (2.6), we have to specify the modal functions $\unicode[STIX]{x1D6FC}_{mn}(t)$ and $\unicode[STIX]{x1D6FD}_{mn}(t)$. Linear model theories show that the modal functions always obey a multidimensional system of ordinary differential equations. Faltinsen & Timokha (Reference Faltinsen and Timokha2009, equation (5.155)) have derived the modal equations in the case of one fluid layer, which are mathematically very similar to our interfacial sloshing problem. By following their procedures, we find the modal equations

(2.8a)$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{\unicode[STIX]{x1D6FC}}_{1n}(t)+2\unicode[STIX]{x1D706}_{1n}\dot{\unicode[STIX]{x1D6FC}}_{1n}(t)+\unicode[STIX]{x1D714}_{1n}^{2}\unicode[STIX]{x1D6FC}_{1n}(t)-\unicode[STIX]{x1D701}_{n}\sin (\unicode[STIX]{x1D6FA}t)=0, & \displaystyle\end{eqnarray}$$

(2.8b)$$\begin{eqnarray}\displaystyle & \displaystyle \ddot{\unicode[STIX]{x1D6FD}}_{1n}(t)+2\unicode[STIX]{x1D706}_{1n}\dot{\unicode[STIX]{x1D6FD}}_{1n}(t)+\unicode[STIX]{x1D714}_{1n}^{2}\unicode[STIX]{x1D6FD}_{1n}(t)+\unicode[STIX]{x1D701}_{n}\cos (\unicode[STIX]{x1D6FA}t)=0, & \displaystyle\end{eqnarray}$$

(2.9)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{1n}^{2}={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}+\unicode[STIX]{x1D6FE}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}\right)^{3}}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}, & \displaystyle\end{eqnarray}$$

(2.10)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}_{n}={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})Rd_{s}\unicode[STIX]{x1D6FA}^{3}}{\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)\right)(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}, & \displaystyle\end{eqnarray}$$

to be fulfilled for $\unicode[STIX]{x1D6FC}_{mn}(t)$ and $\unicode[STIX]{x1D6FD}_{mn}(t)$, where we have included linear damping rates $\unicode[STIX]{x1D706}_{1n}$ expressing the energy dissipation as explained by Faltinsen & Timokha (Reference Faltinsen and Timokha2009, chap. 6). Here, $\unicode[STIX]{x1D714}_{1n}^{2}$ are the natural eigenfrequencies of two-layer gravity–capillary waves in cylindrical tanks, while $\unicode[STIX]{x1D701}_{n}$ can be understood as a mode-dependent forcing parameter describing the strength of excitation. Please note that only non-axisymmetric wave modes $(m=1)$ are excited in the linear regime, possessing exactly one nodal circle. Linear solutions do not exist for $m\neq 1$.

Equations (2.8a) and (2.8b) have the following stationary solutions:

(2.11a)$$\begin{eqnarray}\displaystyle \hspace{-12.0pt}\unicode[STIX]{x1D6FC}_{1n}(t) & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D701}_{n}(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\sin (\unicode[STIX]{x1D6FA}t)-{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D6FA}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\cos (\unicode[STIX]{x1D6FA}t),\qquad\end{eqnarray}$$

(2.11b)$$\begin{eqnarray}\displaystyle \hspace{-12.0pt}\unicode[STIX]{x1D6FD}_{1n}(t) & = & \displaystyle -{\displaystyle \frac{\unicode[STIX]{x1D701}_{n}(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\cos (\unicode[STIX]{x1D6FA}t)-{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D6FA}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\sin (\unicode[STIX]{x1D6FA}t).\qquad\end{eqnarray}$$

(2.12a)$$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-9.60004pt}\unicode[STIX]{x1D719}_{1}(r,\unicode[STIX]{x1D711},z,t)=\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})Rd_{s}\unicode[STIX]{x1D6FA}^{3}}{\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)\right)}}{\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}(z-h_{1})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)}}{\displaystyle \frac{J_{1}\left(\unicode[STIX]{x1D716}_{1n}{\displaystyle \frac{r}{R}}\right)}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}\nonumber\\ \displaystyle & & \displaystyle \hspace{-9.60004pt}\quad \times \left[{\displaystyle \frac{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\sin (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})-{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D6FA}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})\right]\!,\end{eqnarray}$$

(2.12b)$$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-9.60004pt}\unicode[STIX]{x1D719}_{2}(r,\unicode[STIX]{x1D711},z,t)=\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{-(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})Rd_{s}\unicode[STIX]{x1D6FA}^{3}}{\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)\right)}}{\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}(z+h_{2})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}{\displaystyle \frac{J_{1}\left(\unicode[STIX]{x1D716}_{1n}{\displaystyle \frac{r}{R}}\right)}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}\nonumber\\ \displaystyle & & \displaystyle \hspace{-9.60004pt}\quad \times \left[{\displaystyle \frac{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\sin (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})-{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D6FA}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})\right]\!.\end{eqnarray}$$

$\unicode[STIX]{x1D702}$

(2.13)$$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-36.0pt}\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D711},t)=\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})d_{s}R\unicode[STIX]{x1D6FA}^{2}}{\left[(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}\right)^{2}\unicode[STIX]{x1D6FE}\right]}}{\displaystyle \frac{J_{1}\left(\unicode[STIX]{x1D716}_{1n}{\displaystyle \frac{r}{R}}\right)}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}\nonumber\\ \displaystyle & & \displaystyle \hspace{-42.0pt}\quad \times \left[{\displaystyle \frac{\unicode[STIX]{x1D714}_{1n}^{2}(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})+{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D714}_{1n}^{2}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}}\sin (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})\!\right]\!.\end{eqnarray}$$

As a main difference from the previous inviscid theories, the azimuthal part of our solution is composed of both a symmetric ${\sim}\cos (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})$ and an antisymmetric ${\sim}\sin (\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x1D711})$ contribution, whereas inviscid solutions are always symmetric. In the limit $\unicode[STIX]{x1D70C}_{2},\unicode[STIX]{x1D706}_{1n},\unicode[STIX]{x1D6FE}\mapsto 0$ this solution is equivalent to the inviscid one-layer theory by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), as we prove in appendix B. From (2.13) we can deduce the resonance frequencies $\unicode[STIX]{x1D714}_{R,n}$ of the wave modes which are modified by damping. Close to resonance ($\unicode[STIX]{x1D6FA}\approx \unicode[STIX]{x1D714}_{1n}$) the interface elevation is described only by the antisymmetric part of the solution. Therefore

(2.14)$$\begin{eqnarray}\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D711}=\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x03C0}/2,t)|_{\unicode[STIX]{x1D6FA}\approx \unicode[STIX]{x1D714}_{1n}}=A{\displaystyle \frac{2\unicode[STIX]{x1D706}_{1n}\unicode[STIX]{x1D6FA}^{3}\unicode[STIX]{x1D714}_{1n}^{2}}{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D6FA}^{2})^{2}+4\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D6FA}^{2}}},\end{eqnarray}$$

where $A$ contains the frequency-independent prefactors of (2.13). The shaking frequency causing the highest amplitude ($\unicode[STIX]{x1D6FA}\equiv \unicode[STIX]{x1D714}_{R,n}$) is calculated by solving

(2.15)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}(r,\unicode[STIX]{x1D711}=\unicode[STIX]{x1D6FA}t-\unicode[STIX]{x03C0}/2,t)}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}}=0\end{eqnarray}$$

giving a resonance frequency of

(2.16)$$\begin{eqnarray}\unicode[STIX]{x1D714}_{R,n}^{2}=2\sqrt{(\unicode[STIX]{x1D714}_{1n}^{2}-\unicode[STIX]{x1D706}_{1n}^{2})^{2}+\unicode[STIX]{x1D706}_{1n}^{2}\unicode[STIX]{x1D714}_{1n}^{2}}-\unicode[STIX]{x1D714}_{1n}^{2}+2\unicode[STIX]{x1D706}_{1n}^{2}\geqslant \unicode[STIX]{x1D714}_{1n}^{2}.\end{eqnarray}$$

Thus, high damping rates always increase the resonance frequency. This behaviour might appear counterintuitive and is contrary to that of the classical driven harmonic oscillator. It is caused by differences of the external forcing. The forcing amplitude itself depends on the driving frequency and scales with ${\sim}\unicode[STIX]{x1D6FA}^{3}$ (see (2.10)), introducing the tendency that the wave is generally amplified for higher shaking frequencies independently of the resonant amplifications. This effect overcompensates the common frequency drop observed in systems underlying a constant driving force ($\unicode[STIX]{x1D701}_{n}$ independent of $\unicode[STIX]{x1D6FA}$). However, in practical applications, as with OSBs, a noticeable increase of the eigenfrequency will not be observable for the first modes due to their weak dependency on $\unicode[STIX]{x1D706}_{1n}$ in (2.16) for common damping rates $\unicode[STIX]{x1D706}_{1n}\lesssim 1$. An increase of $\unicode[STIX]{x1D714}_{R,n}$ would be observable only in highly overdamped set-ups. For weakly damped systems we identify $\unicode[STIX]{x1D714}_{R,n}\approx \unicode[STIX]{x1D714}_{1n}$ (at least for the first important large-scale wave modes).

2.3 Viscous damping

Equation (2.13) is the final result of the modal analysis containing the responding resonance dynamics of the interface. In contrast to previous inviscid theories (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014; Bouvard et al. Reference Bouvard, Herreman and Moisy2017), our solution does not contain any singularities at the resonance frequencies. Those are resolved by the damping parameters $\unicode[STIX]{x1D706}_{1n}$ determining the maximum amplitudes at resonance $\unicode[STIX]{x1D6FA}\approx \unicode[STIX]{x1D714}_{1n}$. However, using an irrotational approach, the damping rates cannot be further deduced since viscous dissipation is caused by the rotational part of the flow manifested in the boundary layers. Therefore, damping rates must be determined a priori, i.e. by fitting the exponential decay of resonant waves after the shaker is switched off. Then, equation (2.13) can predict the wave amplitudes (within its limits) for all shaking frequencies $\unicode[STIX]{x1D6FA}$ around the considered resonance frequency.

However, to allow for further analysis and a better understanding of how viscous damping can affect the resonance dynamics, it is desirable to describe the wave elevation in direct dependence of the viscosities $\unicode[STIX]{x1D708}_{i}$. For that purpose, we exploit some of the recent results given by Herreman et al. (Reference Herreman, Nore, Guermond, Cappanera, Weber and Horstmann2019), who applied a perturbation approach to study interfacial wave damping. They derived viscous damping rates of free gravity–capillary waves by explicitly calculating Stokes boundary layers for the same two-layer geometry that we are considering in the present paper. It was found in this study that viscous damping rates of free-surface and interfacial waves are inherently different, so that we must carefully distinguish between the one- and two-layer limit of viscous damping rates in the following.

In the first order, the two-layer damping rate $\unicode[STIX]{x1D706}_{1n}^{2L}$ is composed of three different contributions

(2.17)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{1n}^{2L}=\unicode[STIX]{x1D706}_{1n}^{BL}+\unicode[STIX]{x1D706}_{1n}^{IL}+\unicode[STIX]{x1D706}_{1n}^{Int},\end{eqnarray}$$

with

(2.18)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{1n}^{BL} & = & \displaystyle \mathop{\sum }_{i=1,2}\unicode[STIX]{x1D70C}_{i}\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D708}_{i}\unicode[STIX]{x1D6FA}}{8R^{2}}}}\!\left[{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}\!\left(1-{\displaystyle \frac{h_{i}}{R}}\right)\sinh ^{-2}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{i}\right)+\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}^{2}+1}{\unicode[STIX]{x1D716}_{1n}^{2}-1}}\right)\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{i}\right)}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}\!\right]\!\!,\qquad\end{eqnarray}$$

(2.19)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{1n}^{IL} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}}}\sqrt{{\displaystyle \frac{8R^{2}}{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}}}\sqrt{{\displaystyle \frac{8R^{2}}{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D708}_{2}}}}\right)}}\left[{\displaystyle \frac{\left(\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)\right)^{2}}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}\right]\!,\end{eqnarray}$$

(2.20)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{1n}^{Int} & = & \displaystyle {\displaystyle \frac{2\unicode[STIX]{x1D716}_{1n}^{2}\left({\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D708}_{2}}{R^{2}}}-{\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D708}_{1}}{R^{2}}}\right)}{\left({\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\right)}}\left[{\displaystyle \frac{{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{1}\sqrt{\unicode[STIX]{x1D708}_{1}}}}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{2}\sqrt{\unicode[STIX]{x1D708}_{2}}}}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)}{\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)\right)}}\right]\!.\qquad\end{eqnarray}$$

The first contribution $\unicode[STIX]{x1D706}_{1n}^{BL}$ involves viscous dissipation arising at the solid tank boundary layers including the side, top and bottom walls. The second term $\unicode[STIX]{x1D706}_{1n}^{IL}$ describes the dissipation rate in the interfacial boundary layers above and below the interface. Finally, $\unicode[STIX]{x1D706}_{1n}^{Int}$ is the interior damping rate which can be destabilising. Please note that this damping rate is fundamentally different from the well-known damping rates of free-surface waves in upright cylinders calculated by Case & Parkinson (Reference Case and Parkinson1956) and Miles & Henderson (Reference Miles and Henderson1998). The first solid boundary layer term $\unicode[STIX]{x1D706}_{1n}^{BL}$ is physically the same as for one-layer systems. The one-layer limit $\unicode[STIX]{x1D70C}_{1}\mapsto 0$ of $\unicode[STIX]{x1D706}_{1n}^{BL}$ is equivalent to the boundary layer term derived by Case & Parkinson (Reference Case and Parkinson1956). However, the interfacial layer contribution $\unicode[STIX]{x1D706}_{1n}^{IL}$ does not exist in free-surface waves in the leading order. The reason is that the flow field under free surfaces is to a good approximation irrotational, which is not the case around moving interfaces, where viscous boundary layers evolve to balance the strong tangential shear flows between the liquids. Indeed, we find $\unicode[STIX]{x1D706}_{1n}^{IL}\mapsto 0$ for $\unicode[STIX]{x1D70C}_{1}\mapsto 0$. For liquids with similar densities and viscosities $\unicode[STIX]{x1D706}_{1n}^{BL}$ and $\unicode[STIX]{x1D706}_{1n}^{IL}$ can have comparable magnitudes and should always both be considered. In contrast, the interior damping rate $\unicode[STIX]{x1D706}_{1n}^{Int}$ is completely negligible for liquids with densities of the same order. However, this term cannot be ignored anymore when we approach the free-surface limit $\unicode[STIX]{x1D70C}_{1}\mapsto 0$. Then, the interior damping is increased by orders of magnitude until it reaches the limit $\unicode[STIX]{x1D706}_{1n}^{Int}\mapsto 2\unicode[STIX]{x1D708}_{2}\unicode[STIX]{x1D716}_{1n}^{2}/R^{2}$, the familiar interior damping rate of free-surface waves that is generally not negligible. This is the only remaining dissipation term for irrotational surface waves, while interfacial wave damping is always manifested in rotational boundary layers. All in all, leading-order interfacial damping is well described by

(2.21)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{1n}^{2L}=\unicode[STIX]{x1D706}_{1n}^{BL}+\unicode[STIX]{x1D706}_{1n}^{IL}.\end{eqnarray}$$

In contrast, free-surface damping must be calculated using the contributions

(2.22)$$\begin{eqnarray}\unicode[STIX]{x1D706}_{1n}^{1L}=\unicode[STIX]{x1D706}_{1n}^{BL}+\tilde{\unicode[STIX]{x1D706}}_{1n}^{Int},\end{eqnarray}$$

with

(2.23)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}_{1n}^{BL}=\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D708}_{2}\unicode[STIX]{x1D6FA}}{8R^{2}}}}\left[{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}^{2}+1}{\unicode[STIX]{x1D716}_{1n}^{2}-1}}+{\displaystyle \frac{2\unicode[STIX]{x1D716}_{1n}\left(1-{\displaystyle \frac{h_{2}}{R}}\right)}{\sinh \left(2{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}\right]\!, & \displaystyle\end{eqnarray}$$

(2.24)$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D706}}_{1n}^{Int}={\displaystyle \frac{\unicode[STIX]{x1D708}_{2}}{R^{2}}}\left[2\unicode[STIX]{x1D716}_{1n}^{2}-{\displaystyle \frac{1}{\unicode[STIX]{x1D716}_{1n}^{2}-1}}\left(1+{\displaystyle \frac{2\unicode[STIX]{x1D716}_{1n}h_{2}}{R\sinh \left(2{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}}{R}}h_{2}\right)}}\right)\right]\!. & \displaystyle\end{eqnarray}$$

The boundary layer damping (2.23) is the one-layer limit of (2.18), whereas the interior damping rate $\tilde{\unicode[STIX]{x1D706}}_{1n}^{Int}$ was supplemented by an additional sidewall contribution calculated by Miles & Henderson (Reference Miles and Henderson1998) in order to describe damping with the best known accuracy. This contribution is entirely negligible for interfacial waves and still small for typical free-surface cylinders such as OSBs. In this form, the damping rate (2.22) is fully equivalent to the formulation recently used by Raynovskyy & Timokha (Reference Raynovskyy and Timokha2018a) to study resonant sloshing. In the shallow water limit, the damping rate is further consistent with the Stokes wall shear stress model presented by Alpresa et al. (Reference Alpresa, Sherwin, Weinberg and van Reeuwijk2018b), predicting a spiral-like shape of the horizontal velocity profile in the bottom boundary layer.

Technically, both damping rates $\unicode[STIX]{x1D706}_{1n}^{2L}$ and $\unicode[STIX]{x1D706}_{1n}^{1L}$ are only valid for free gravity–capillary waves; dissipation rates of forced wave motion are generally more complex. However, damping affects linear waves mainly in a small window around the eigenfrequencies $\unicode[STIX]{x1D6FA}\approx \unicode[STIX]{x1D714}_{1n}$, as we will show in § 2.4. In this frequency range, the interfacial motions are very similar to free wave motions and dissipation rates can be well approximated by (2.17) and (2.22). This is not generally the case for frequencies far below or above the resonance frequency, but there, damping is of no importance except for a largely overdamped system.

Finally, it is important to note that the damping rates do not involve any dissipation mechanism associated with the contact line (contact line hysteresis, meniscus dynamics) or surface contamination. Hence, (2.17) and (2.22) are always expected to slightly underestimate viscous damping rates and should be considered as conservative estimations to predict maximum expectable interface elevations. While it is often difficult to realise two-layer stratifications involving free-sliding contact lines and negligible menisci (Horstmann et al. Reference Horstmann, Wylega and Weier2019), these idealised boundary conditions apply for most mid- and large-size free-surface systems, including typical commercial OSBs.

2.4 Dimensionless prediction of resonance dynamics

In this section we want to relate the wave solutions (2.12b) and (2.13) to the eight dimensionless parameters introduced in § 2.1 in order to discuss the damped wave dynamics predicted by our model. We can write the solutions in dimensionless forms by introducing the dimensionless variables $\tilde{r}=r/R,~\tilde{z}=z/R$ and $\tilde{t}=\unicode[STIX]{x1D6FA}t$. This way, the interface elevation $\unicode[STIX]{x1D702}$ (2.13) can be expressed as

(2.25)$$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D702}(\tilde{r},\unicode[STIX]{x1D711},\tilde{t})}{R}}=\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{2Fr}{\left(1+{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}^{2}}{Bo}}\right)}}{\displaystyle \frac{J_{1}(\unicode[STIX]{x1D716}_{1n}\tilde{r})}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \left[{\displaystyle \frac{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}}{{\displaystyle \frac{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)^{2}}{\unicode[STIX]{x1D6EC}_{1n}}}+\unicode[STIX]{x1D6EC}_{1n}}}\sin (\,\tilde{t}-\unicode[STIX]{x1D711})+{\displaystyle \frac{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)+{\displaystyle \frac{\unicode[STIX]{x1D6EC}_{1n}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)}}}}\cos (\,\tilde{t}-\unicode[STIX]{x1D711})\right]\!,\qquad\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}(Fr,E,Bo,H_{i},A)=\unicode[STIX]{x1D714}_{1n}/\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D6EC}_{1n}(Re_{i},H_{i},A)=2\unicode[STIX]{x1D706}_{1n}/\unicode[STIX]{x1D6FA}$. Full formulas for these dimensionless frequencies and damping rates as well as for the dimensionless potentials and velocities are given in appendix A. Equation (2.25) is the final result of our theoretical work, allowing analysis of the wave elevation as a function of all eight dimensionless input parameters. The physics contained within this solution is very rich, although it may be difficult to recognise in the presented form. By analysing the different wave contributions and the asymptotics, we can further deduce various wave properties. One main difference between this damped wave solution in comparison with the previous inviscid description (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) is the existence of an antisymmetric part ${\sim}\sin (\,\tilde{t}-\unicode[STIX]{x1D711})$, not present in the inviscid solution. The superposition of the symmetric and antisymmetric solutions captures the phase shift $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ between the shaking table and the responding wave. The phase $\unicode[STIX]{x1D711}^{\ast }$ of the maximum wave elevation can be calculated by solving

(2.26)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}(\tilde{r},\unicode[STIX]{x1D711}^{\ast },\tilde{t})}{R\unicode[STIX]{x2202}\unicode[STIX]{x1D711}^{\ast }}}=0.\end{eqnarray}$$

From this, we find a phase difference of

(2.27)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(\tilde{r})=\text{arctan}\left({\displaystyle \frac{\displaystyle \mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{1}{\left(Bo+\unicode[STIX]{x1D716}_{1n}^{2}\right)}}{\displaystyle \frac{J_{1}(\unicode[STIX]{x1D716}_{1n}\tilde{r})}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}{\displaystyle \frac{\unicode[STIX]{x1D6EC}_{1n}\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)^{2}+\unicode[STIX]{x1D6EC}_{1n}^{2}}}}{\displaystyle \mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{1}{\left(Bo+\unicode[STIX]{x1D716}_{1n}^{2}\right)}}{\displaystyle \frac{J_{1}(\unicode[STIX]{x1D716}_{1n}\tilde{r})}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}{\displaystyle \frac{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)}{\left(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1\right)^{2}+\unicode[STIX]{x1D6EC}_{1n}^{2}}}}}\right).\end{eqnarray}$$

The phase difference depends on the radial coordinate $\tilde{r}$ since the crests of the individual nodal cycles $n$, which are located at different radial positions, can be shifted differently by the mode-dependent damping rates (2.17). Considering only the phase shift of the first mode ($n=1$), equation (2.27) simplifies to

(2.28)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}|_{n=1}=\text{arctan}\left({\displaystyle \frac{\unicode[STIX]{x1D6EC}_{11}}{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1}}\right)=\text{arctan}\left({\displaystyle \frac{2\unicode[STIX]{x1D706}_{11}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x1D714}_{11}^{2}-\unicode[STIX]{x1D6FA}^{2}}}\right),\end{eqnarray}$$

exactly reproducing the well-known phase shift of the driven harmonic oscillator, even though we found a different resonance frequency (see (2.16)). Equation (2.28) can be exploited to estimate the linear out-of-phase flow. By defining an expedient critical small phase shift $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}>\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{c}$ as an out-of-phase threshold, a corresponding critical shaking frequency $\unicode[STIX]{x1D6FA}_{c}$ for out-of-phase flow can be calculated as

(2.29)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D6FA}_{c}}{\unicode[STIX]{x1D714}_{11}}}\geqslant \sqrt{{\displaystyle \frac{1}{\unicode[STIX]{x1D6EC}_{11}\cot (\unicode[STIX]{x0394}\unicode[STIX]{x1D711}_{c})+1}}}.\end{eqnarray}$$

It becomes apparent that the linear phase shift is primarily specified by the damping forces.

Likewise, the maximum achievable interface elevation for given system parameters can be estimated from (2.25). Close to resonance the wave elevation is determined only by the antisymmetric part of the solution (2.25). Here, we can locate the maximum elevation $\unicode[STIX]{x1D702}_{0}/R$ of the first dominant mode at $\tilde{r}=1$, $\unicode[STIX]{x1D711}=\tilde{t}-\unicode[STIX]{x03C0}/2$ and $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D714}_{11}$, yielding

(2.30)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}={\displaystyle \frac{2Fr}{\left(1+{\displaystyle \frac{\unicode[STIX]{x1D716}_{1n}^{2}}{Bo}}\right)}}{\displaystyle \frac{1}{\unicode[STIX]{x1D716}_{11}^{2}-1}}{\displaystyle \frac{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}}{\unicode[STIX]{x1D6EC}_{11}}}.\end{eqnarray}$$

This expression is only valid for finite amplitudes. Nonlinear effects and wave breaking cause additional dissipation, effectively reducing the maximal elevation, an effect not captured by the theory. However, equation (2.30) is a good measure for the maximum attainable wave elevation. As pointed out by Reclari (Reference Reclari2013) and Alpresa et al. (Reference Alpresa, Sherwin, Weinberg and van Reeuwijk2018a), linear potential theory can adequately predict free-surface motions up to a critical amplitude, where breaking is expected. This usually happens if the wave steepness exceeds a critical value beyond which gravity waves become unstable. For the first and practically most important wave mode $\unicode[STIX]{x1D716}_{11}$, the critical wave steepness $\unicode[STIX]{x1D702}_{0}/R\gtrsim 0.44$ has been well established as a wave breaking criterion in orbitally shaken cylinders with deep fluid layers (Reclari Reference Reclari2013). Up to this wave slope, equation (2.30) is expected to predict the resonant interfacial displacement with sufficient accuracy. However, this criterion is not valid for shallow layers, where the vicinity of the bottom leads to additional effects including a partially uncovered bottom. For these cases Alpresa et al. (Reference Alpresa, Sherwin, Weinberg and van Reeuwijk2018a) verified the shallow breaking criterion $\unicode[STIX]{x1D702}_{0}/R\gtrsim 0.7H$. By combining these criteria as upper thresholds and by modifying (2.30) according to (A 4) and (A 5), we can derive the following relations for the resonant wave slope within the one-layer limit as relevant for OSBs:

(2.31)$$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}=\unicode[STIX]{x1D712}\sqrt{Re}E\text{ if }\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}<0.44 & \text{for }H\gtrsim 0.63,\\[7.0pt] {\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}<0.7H & \text{for }H\lesssim 0.63,\end{array}\right.\nonumber\\ \displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}\leqslant \unicode[STIX]{x1D712}\sqrt{Re}E\text{ if }\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}\geqslant 0.44 & \text{for }H\gtrsim 0.63,\\[7.0pt] {\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}\geqslant 0.7H & \text{for }H\lesssim 0.63,\end{array}\right.\end{eqnarray}$$

with the form factors

$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D712}={\displaystyle \frac{\sqrt{8}\unicode[STIX]{x1D716}_{11}\sinh ^{2}(\unicode[STIX]{x1D716}_{11}H)}{\unicode[STIX]{x1D716}_{11}(\unicode[STIX]{x1D716}_{11}^{2}-1)(1-H)+\sinh (\unicode[STIX]{x1D716}_{11}H)\cosh (\unicode[STIX]{x1D716}_{11}H)(\unicode[STIX]{x1D716}_{11}^{2}+1)}}.\nonumber\end{eqnarray}$$

Here, the index $i=2$ was omitted for better readability. In line with a conservative estimation, the interior damping distribution ${\sim}Re$ in (2.22) was ignored since it is, in comparison with boundary layer damping, negligibly small for lower wave modes in OSB-size cylinders ($R\sim 10~\text{cm}$). The maximum wave displacement is independent of the gravity force and surface tension. It is solely determined by the balance of driving and viscous forces. We can determine the relation $\sqrt{Re}E=d_{s}/(2\unicode[STIX]{x1D6FF})$, showing that the maximum slope scales approximately with the length ratio between the shaking radius $d_{s}/2$ and the typical boundary layer thickness $\unicode[STIX]{x1D6FF}$. This scaling is precisely valid for interfacial waves, where interior damping is always negligible. The implied scaling law $\unicode[STIX]{x1D702}_{0}/R\sim Re^{\unicode[STIX]{x1D705}}$ with $\unicode[STIX]{x1D705}=0.5$ can generally be considered as an upper threshold also for non-resonant excitations since, at resonance, the impact of damping forces is always at a maximum. We finally conclude that an upper value of the scaling coefficient of $\unicode[STIX]{x1D705}\leqslant 0.5$ can be expected for most OSBs and all two-layer systems depending on the shaking frequency $\unicode[STIX]{x1D6FA}$ and on nonlinearities.

Lastly, we want to analyse the wave character of the elevated free surface. For low shaking frequencies $\unicode[STIX]{x1D6FA}$ sufficiently below the first eigenfrequency $\unicode[STIX]{x1D714}_{11}$ the displaced free surface has the form of a tilted disk rather than a wavy shape (Reclari Reference Reclari2013; Weheliye et al. Reference Weheliye, Yianneskis and Ducci2013; Ducci & Weheliye Reference Ducci and Weheliye2014). The transition from the flat disk shape to the waveform is important for OSBs since the mixing dynamics is changed considerably: as soon as the surface becomes wavy, the nonlinear out-of-phase flow regime starts to develop (Weheliye et al. Reference Weheliye, Yianneskis and Ducci2013).

As shown in appendix B, the symmetric part of solution (2.25) can be transformed into the following form:

(2.32)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D702}(\tilde{r},\unicode[STIX]{x1D711},\tilde{t})}{R}}=2Fr\left[{\displaystyle \frac{\tilde{r}}{2}}+\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)-\unicode[STIX]{x1D6EC}_{1n}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{1n}}^{2}-1)^{2}+\unicode[STIX]{x1D6EC}_{1n}^{2}}}{\displaystyle \frac{J_{1}(\unicode[STIX]{x1D716}_{1n}\tilde{r})}{(\unicode[STIX]{x1D716}_{1n}^{2}-1)J_{1}(\unicode[STIX]{x1D716}_{1n})}}\right]\cos (\,\tilde{t}-\unicode[STIX]{x1D711}),\end{eqnarray}$$

when surface tension is neglected ($Bo\gg 1$). The radial part of (2.32) is composed of a linear contribution ${\sim}\tilde{r}$ representing the disk and a wavy contribution ${\sim}J_{1}(\unicode[STIX]{x1D716}_{1n}\tilde{r})$ resulting from the superposition of free gravity wave modes. We can now compare the linear contribution with the wave-like contributions. We consider a fixed time point $\tilde{t}=0$ and the radial position $\tilde{r}=1$, where both the disc and the first wave mode are of maximum elevation. Further restricting ourselves to small shaking frequencies $\unicode[STIX]{x1D6FA}$ below the first resonance frequency $\unicode[STIX]{x1D714}_{11}$ allows us to keep only the first wave mode ($n=1$). We can then compare the wavy parts (mode $n=1$ in (2.32) and the antisymmetric part of (2.25)) with the disc component and calculate the point where the sum of the wavy components matches a certain percentage $\unicode[STIX]{x1D6F6}$ of the disc solution

(2.33)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D6F6}}{2}}={\displaystyle \frac{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1)-\unicode[STIX]{x1D6EC}_{11}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1)^{2}+\unicode[STIX]{x1D6EC}_{11}^{2}}}{\displaystyle \frac{1}{(\unicode[STIX]{x1D716}_{11}^{2}-1)}}-{\displaystyle \frac{\unicode[STIX]{x1D6EC}_{11}\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1)^{2}+\unicode[STIX]{x1D6EC}_{11}^{2}}}{\displaystyle \frac{1}{(\unicode[STIX]{x1D716}_{11}^{2}-1)}}\tan (\unicode[STIX]{x1D711}^{\ast }).\end{eqnarray}$$

The phase $\unicode[STIX]{x1D711}^{\ast }$ of maximum wave elevation can be calculated in analogy to (2.27) and is given by

(2.34)$$\begin{eqnarray}\unicode[STIX]{x1D711}^{\ast }=\arctan \left(-{\displaystyle \frac{\unicode[STIX]{x1D6EC}_{11}}{\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1}}\right).\end{eqnarray}$$

By inserting (2.34) into (2.33) the wave percentage simplifies to

(2.35)$$\begin{eqnarray}\unicode[STIX]{x1D6F6}={\displaystyle \frac{2}{(\unicode[STIX]{x1D716}_{11}^{2}-1)}}{\displaystyle \frac{1}{(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D714}_{11}}^{2}-1)}}.\end{eqnarray}$$

The damping parameter $\unicode[STIX]{x1D6EC}_{11}$ does not occur in (2.35) and thus the wave fraction $\unicode[STIX]{x1D6F6}$ is independent of the system’s damping. This somewhat surprising result implies that the disc and the first wave mode are equally damped. Hence, damping will not influence the plane-wave transition. Equation (2.35) can be derived as well by directly using the inviscid solution (equation (2.32) with $\unicode[STIX]{x1D6EC}_{11}=0$). Any possible influence of viscosity to the transition would be manifested in the rotational part of the flow confined to the boundary layers.

By including the dimensionless frequency (A 4) in (2.35) we can derive the critical Froude number $Fr_{wave}$ needed to induce a certain waviness $\unicode[STIX]{x1D6F6}$

(2.36)$$\begin{eqnarray}Fr_{wave}(\unicode[STIX]{x1D6F6})={\displaystyle \frac{E\unicode[STIX]{x1D716}_{11}\tanh (\unicode[STIX]{x1D716}_{11}H)}{2(\unicode[STIX]{x1D716}_{11}^{2}-1)^{-1}\unicode[STIX]{x1D6F6}^{-1}+1}}.\end{eqnarray}$$

This formula can be used to predict the limits of the linear scaling law $\unicode[STIX]{x1D702}_{0}/R\sim Fr$ proposed by Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) and Ducci & Weheliye (Reference Ducci and Weheliye2014) as well as the nonlinear out-of-phase flow regime in OSBs, as we will discuss in § 4.3.

2.5 Limits of applicability

Two of the initial assumptions are mainly responsible for the limits of our theoretical model: irrotationality and linearity. Even though the dissipative effect of rotational boundary layers is incorporated in our model by means of the damping rate formula (2.17), we neglect the influences of the boundary layers on the wave shapes and the velocity fields in the bulk since the flow solutions violate the no-slip boundary conditions. Hence, the wave modes (2.25) can describe the bulk flow only adequately if the boundary layer thicknesses $\unicode[STIX]{x1D6FF}_{i}$ are far smaller than the internal dimensions

(2.37)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{i}\sim \sqrt{{\displaystyle \frac{\unicode[STIX]{x1D708}_{i}}{\unicode[STIX]{x1D6FA}}}}={\displaystyle \frac{R}{\sqrt{Re_{i}}}}\ll R.\end{eqnarray}$$

This estimate directly leads to the condition

(2.38)$$\begin{eqnarray}\sqrt{Re_{i}}\gg 1,\end{eqnarray}$$

which is no serious limitation in practice. Even small-scale bioreactors rarely operate at $Re\lesssim 10^{2}$ (Alpresa et al. Reference Alpresa, Sherwin, Weinberg and van Reeuwijk2018a) since the fluid-driving wave would otherwise be too small to generate sufficient mixing.

A more serious restriction is the linearisation limiting our model to fairly small amplitude waves. Nonlinear effects such as subharmonic wave modes (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), the mean flow (Bouvard et al. Reference Bouvard, Herreman and Moisy2017) and turbulent wave dynamics are not captured. As discussed before, the wave breaking conditions

(2.39)$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}\lesssim \left\{\begin{array}{@{}ll@{}}0.44, & \text{for }\min (H_{1},H_{2})\gtrsim 0.63,\\[2.0pt] 0.7\min (H_{1},H_{2}), & \text{for }\min (H_{1},H_{2})\lesssim 0.63,\end{array}\right. & & \displaystyle\end{eqnarray}$$

are good estimates for linearity since good agreement with the inviscid potential model was found for the single-crested $n=1$ mode in this amplitude regime. For OSBs we can use (2.31) to derive an upper threshold for the Reynolds number

(2.40)$$\begin{eqnarray}\displaystyle \sqrt{Re}\lesssim \left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{0.44}{\unicode[STIX]{x1D712}E}}, & \text{for }H\gtrsim 0.63,\\[12.0pt] {\displaystyle \frac{0.7H}{\unicode[STIX]{x1D712}E}}, & \text{for }H\lesssim 0.63.\end{array}\right. & & \displaystyle\end{eqnarray}$$

Consequently, a critical eccentricity $E$ can be estimated beyond which nonlinear waves might evolve for excitation frequencies lower than or equal to the resonance frequency. In practice, condition (2.40) is fulfilled for small-scale OSBs of order $R\sim 1~\text{cm}$ as used for protein degradation. But it is violated for large-scale bioreactors of orders $R\sim 10$ up to $R\sim 100~\text{cm}$, as commonly used, e.g. for mammalian cell cultivation (Klöckner et al. Reference Klöckner, Tissot, Wurm and Büchs2012, Reference Klöckner, Diederichs and Büchs2014). However, equation (2.40) is a conservative limit for the maximum attainable amplitudes at resonance. For excitation frequencies sufficiently below or above the resonance frequency, amplitudes are considerably lower and our theory can predict wave amplitudes up to a critical Froude number, as previously shown for the inviscid model (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). Some of the results, such as the prediction of the waviness (2.36), are of practical relevance for large-scale OSBs as well; see § 4.3 for details.

An inequality similar to (2.39) applies for the two-layer case, albeit with different wave breaking limits, which can depend on the density difference. Even so, an application to two-layer systems extends the valid parameter range of our model in comparison with equivalent one-layer waves (Horstmann et al. Reference Horstmann, Wylega and Weier2019). Here, the model performs best. This is mainly due to the density differences, which are considerably lower for most liquid–liquid interfaces compared with gas–liquid free-surfaces. Typical oil–water combinations possess up to an order of magnitude lower eigenfrequencies compared with free surfaces. In consequence, the wave amplitudes for frequencies near the resonance frequency, where the elevation scales with ${\sim}\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D714}_{1n}^{2}$ (cf. (2.13)), are much smaller.

Concluding this section, we have to emphasise again that our theory cannot describe contact line dynamics. The theory is valid only for idealised boundary conditions (free-sliding contact line with a static contact angle of $90^{\circ }$). This condition is in good approximation fulfilled for many mid- and large-scale surface-wave systems but harder to realise for interfacial waves (cf. Horstmann et al. Reference Horstmann, Wylega and Weier2019). For the latter case, contact line hysteresis and dynamic menisci can frequently be observed and pose great challenges to modelling. For the idealised case of a perfectly fixed contact line there are some wave theories in the literature (Miles Reference Miles1991; Henderson & Miles Reference Henderson and Miles1994). Very recently, Viola & Gallaire (Reference Viola and Gallaire2018) succeeded in developing a mathematical framework which incorporates a dynamic contact angle model. Formulating these approaches under orbital excitation is a promising task that we have to leave for future studies.

4.1 Two-layer wave elevation comparisons with Horstmann et al. (Reference Horstmann, Wylega and Weier2019)

The multi-layer orbital sloshing experiment by Horstmann et al. (Reference Horstmann, Wylega and Weier2019) was designed with the intention to imitate the wave dynamics excited by the metal pad roll instability in liquid metal batteries (Weber et al. Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017; Horstmann et al. Reference Horstmann, Weber and Weier2018). Controlling the contact line dynamics was the major experimental difficulty. For most liquid combinations Horstmann et al. (Reference Horstmann, Wylega and Weier2019) observed entirely or partially fixed contact lines, with the latter additionally subject to a contact angle hysteresis. Complex contact line dynamics affects both the wave modes and the damping rates and is not covered in our model, however, recently, Viola & Gallaire (Reference Viola and Gallaire2018) developed a theoretical framework treating all these effects. However, the particular combination of Paraffinum Perliquidum (PP) ($\unicode[STIX]{x1D70C}_{1}=846~\text{g}~\text{cm}^{-3}$, $\unicode[STIX]{x1D708}_{1}=36~\text{mm}^{2}~\text{s}^{-1}$) overlaying Wacker^® AK 35 silicone oil ($\unicode[STIX]{x1D70C}_{2}=955~\text{g}~\text{cm}^{-3}$, $\unicode[STIX]{x1D708}_{2}=35~\text{mm}^{2}~\text{s}^{-1}$) fulfilled the idealised boundary condition we have used in this study: no meniscus was visible and the contact line slid almost freely along the sidewalls. Hence, these oils are a promising choice for comparisons.

Figure 5.

(a) Experimental resonance curves (circular markers) of paraffin oil$|$silicon oil interface elevations at the position $\tilde{r}=0.84$ presented in Horstmann et al. (Reference Horstmann, Wylega and Weier2019) for four chosen values of $H_{2}$ in comparison with the theoretical prediction (2.25). Resonance curves are depicted by the dimensionless local wave amplitudes as a function of the shaking frequency $f=\unicode[STIX]{x1D6FA}/2\unicode[STIX]{x03C0}$. (b) Phase shift $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ between wave and shaker for the same measurements. The theoretical phase shifts were calculated using (2.27).

Horstmann et al. (Reference Horstmann, Wylega and Weier2019) observed that the measured resonance frequencies differed from the theoretical natural eigenfrequencies (2.9) in many stratifications. The inviscid resonance curves appeared shifted relative to the measured ones. The reason for these deviations were, on the one hand, the fixed contact lines in oil$|$water stratifications that increased the peak frequency with respect to (2.9). On the other hand, a similar shift was also observable in many oil$|$oil combinations as for PP$|$AK 35, caused by partial mixing at the interface. Partial mixing led to a temporally increasing reduction of the effective density difference $\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}$, eventually lowering the measured peak frequency with respect to (2.9).

We can nevertheless compare the theory with the experiments when using the measured resonance frequency in (2.25) instead of the theoretical frequency (2.9). In figure 5(a), resonance curves calculated accordingly are compared with the PP$|$AK 35 measurements of Horstmann et al. (Reference Horstmann, Wylega and Weier2019) for four representatively selected values of $H_{2}$. A small cylinder of radius $R=5~\text{cm}$ and a shaking diameter of $d_{s}=2.5~\text{cm}$ were used, and the interface elevation was measured at a fixed position $r=4.2~\text{cm}$. The theoretical prediction coincides well with the experimental data around the first resonance $f\lesssim 60~\text{r.p.m.}$; only the amplitudes close to the peaks are slightly overestimated for the shallow cases. This was expected since additional dissipation mechanisms (e.g. based on interface contamination or the contact line dynamics) are not covered by the theoretical damping rate (2.17). They are, however, unavoidable in the experiments. The direct comparison between measurements and theory in Herreman et al. (Reference Herreman, Nore, Guermond, Cappanera, Weber and Horstmann2019) demonstrates that the theory slightly underpredicts the measured damping rates. Nevertheless, the first peak is fairly well captured even for the highly damped case $H_{2}=0.3$, where, for $f\gtrsim 38~\text{r.p.m.}$, the existing inviscid theory largely fails. The second peak is completely damped out in our experiments. This behaviour is reflected in the theory, although there is a disagreement in the range $60~\text{r.p.m.}\lesssim f\lesssim 80~\text{r.p.m.}$ We observed in the experiment that high excitation frequencies $f\gtrsim 60~\text{r.p.m.}$ lead to interface ripples, probably causing substantial increase of the system’s dissipation and explaining the mismatch here. Qualitatively, the theory can describe how amplitude saturations can result from the superposition of higher overdamped wave modes.

Complementarily, figure 5(b) shows the theoretical and experimental phase shifts $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$ for the same measurements. As for the amplitudes, the phase shifts around the first resonance are well described by our model, while the experimentally observed phase jumps for $H_{2}=0.4$ and $H_{2}=0.3$ are somewhat smoother than predicted due to the slightly underestimated damping discussed before. For higher frequencies, the theoretical phase shifts overshoot the measured curves. However, here, the theory (2.27) qualitatively reveals why a phase shift of $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}=180\,^{\circ }$ is never reached. It predicts individual phase shifts for every eigenfrequency which are progressively smoothed out by increased damping. If the damping is strong enough, as in the present case, the second mode causes a reversed phase shift before a full out-of-phase flow can be developed. Eventually, the superposition of strongly smoothed phase transitions at higher modes yields the nearly linear growth we observe in 5(b) at the highest frequencies. Hence, our comparisons show that the model captures all the relevant physics causing this damped resonance dynamics. Because it is challenging to control the experimental conditions, the experimental damping rates are always slightly underestimated by the theory. However, all measured resonance and phase curves shown here can be almost perfectly reproduced by the theory, if fitted values of the damping rates instead of analytical values (2.21) are used. Following this procedure, our model can also be employed to determine linear damping rates by measuring only wave amplitudes or phase shifts. This is particularly useful when damping is high and the exponential decay happens too rapidly for a precise calculation of the decay rates.

4.2 Velocity comparisons with Bouvard et al. (Reference Bouvard, Herreman and Moisy2017)

Bouvard et al. (Reference Bouvard, Herreman and Moisy2017) created an orbital sloshing experiment with the aim to investigate the wave-induced mean mass transport in the weakly nonlinear regime. They conducted a series of measurements in a glass cylinder of radius $R=5.12~\text{cm}$ and height $h_{2}=11.1~\text{cm}$. Silicon oils with a kinematic viscosity of either $\unicode[STIX]{x1D708}=50$ or $500~\text{mm}^{2}~\text{s}^{-1}$ and surface tension of $\unicode[STIX]{x1D6FE}=21\times 10^{-3}~\text{N}~\text{m}^{-1}$ were used as working fluids. Velocity fields were measured by particle imaging velocimetry at different vertical and horizontal light sheet positions. While the study focuses on the mean flow, which is not covered in our approach, other fundamental flow properties such as the wave fields and phase shifts are reported as well. Bouvard et al. (Reference Bouvard, Herreman and Moisy2017) defined the norm of the horizontal velocity at the cell centre $r=0$ and height $z=z_{0}$ as a measure of the wave amplitude

(4.1)$$\begin{eqnarray}|\mathbf{u}_{\bot }|(r=0,z=z_{0})=\sqrt{u_{r}^{2}+u_{\unicode[STIX]{x1D711}}^{2}}.\end{eqnarray}$$

Figure 6(a) shows the normalised wave amplitude as a function of the normalised shaking frequency $\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D714}_{11}$ at $z_{0}=-0.23R$ under excitation with $E=0.057$ for both silicon oils. The theoretical curves were calculated by the potential solution (A 1b) together with (2.5) in the one-layer limit $\unicode[STIX]{x1D70C}_{1}=0$. It becomes clear that our approach does not provide an advantage for $\unicode[STIX]{x1D708}=50~\text{mm}^{2}~\text{s}^{-1}$. The damped solution mainly coincides with the inviscid theory and highly overestimates the peak velocities. This is not unexpected since Bouvard et al. (Reference Bouvard, Herreman and Moisy2017) report a hysteresis of the resonant wave dynamics and thus clearly nonlinear behaviour. Here, we exceed the limits of our approach. The nonlinear multimodal theory by Faltinsen, Lukovsky & Timokha (Reference Faltinsen, Lukovsky and Timokha2016) was recently adapted to orbital sloshing (Raynovskyy & Timokha Reference Raynovskyy and Timokha2018b) and allows the first predictions of steady-state sloshing dynamics in this resonant regime. An advantage of our model is visible for $\unicode[STIX]{x1D708}=500~\text{mm}^{2}~\text{s}^{-1}$. Here, our theory tends to overestimate the resonant velocities as well. This might be due to the influence of ingrowing boundary layers or additional experimental damping mechanisms not captured by our theory (surface contamination, slight contact line hysteresis). However, for frequencies higher than the resonance frequency, the measured wave velocities agree better with our theory than with the inviscid model.

Figure 6.

(a) Normalised velocity amplitude measured at $r=0$, $z_{0}=-0.23R$ as a function of the dimensionless shaking frequency $\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D714}_{11}$ for both silicon oils and $E=0.057$. The markers show the measurements conducted by Bouvard et al. (Reference Bouvard, Herreman and Moisy2017), the black line our theoretical prediction due to the potential solution (2.12b) and the dashed line the inviscid prediction by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). (b) Phase shifts as a function of $\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D714}_{11}$ for the same measurements. (c) Normalised velocity amplitude as a function of $E$ of the viscous silicon oil ($\unicode[STIX]{x1D708}_{2}=500~\text{mm}^{2}~\text{s}^{-1}$) for three different shaking frequencies.

Generally, it is known that potential models can more accurately predict interface elevations (especially for higher amplitudes) than velocity fields. The latter are influenced to a larger degree by boundary layers and secondary flow structures not captured in potential approaches (Ibrahim Reference Ibrahim2005). This observation is supported by figure 6(b), showing the corresponding phase shifts for both oils. Our model (2.28) can predict the frequency-dependent shift well even for $\unicode[STIX]{x1D708}=50~\text{mm}^{2}~\text{s}^{-1}$, although it failed to predict the resonant velocities. The phase shift of the centre velocity coincides with the phase shift of the surface elevation, which seems to be much more robust towards effects of (weakly) nonlinear waves and secondary rotational flows.

Bouvard et al.’s (Reference Bouvard, Herreman and Moisy2017) measurements of the horizontal centre velocity obtained with the viscous silicon oil $\unicode[STIX]{x1D708}=50~\text{mm}^{2}~\text{s}^{-1}$ for varying shaking diameters $E$ at three different shaking frequencies are shown in figure 6(c). A linear scaling between the velocity and $E$ becomes apparent. For the two measurement series $\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D714}_{11}=0.67$ and 0.78 conducted at frequencies sufficiently below the resonance, Bouvard et al. (Reference Bouvard, Herreman and Moisy2017) found a good agreement with the inviscid theory. The damped predictions scarcely deviate from the inviscid ones and the effect of damping is weak as discussed before. But for the measurement $\unicode[STIX]{x1D6FA}/\unicode[STIX]{x1D714}_{11}=0.89$, which is inside the resonant regime, our damped model coincides considerably better. The curve fits the first three or four measurements, until higher forcing parameters $E$ excite a nonlinear dynamics, manifested in an incipient saturation.

4.3 Scaling law and nonlinear out-of-phase flow prediction in comparison with Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013)

Recently, Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013), Ducci & Weheliye (Reference Ducci and Weheliye2014) and Weheliye et al. (Reference Weheliye, Cagney, Rodriguez, Micheletti and Ducci2018) have made significant progress in identifying and understanding the fluid mechanics of OSBs with a focus on the mixing dynamics. Most OSBs operate in strongly nonlinear regimes, which precludes comparison of our theoretical findings. However, these works reported two different physical effects of high practical relevance, which can be described by our theory under certain assumptions discussed and verified below.

A main result of Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) is the existence of a linear scaling law

(4.2)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D702}_{0}}{R}}=\unicode[STIX]{x1D6FC}_{0}Fr,\end{eqnarray}$$

very useful for estimating surface elevations in OSBs. The constant $\unicode[STIX]{x1D6FC}_{0}$ is reported to depend on the working fluid. Its value is $\unicode[STIX]{x1D6FC}_{0}=1.4$ for water, while lower values were measured for silicon oils of higher viscosity (Ducci & Weheliye Reference Ducci and Weheliye2014). The scaling law was experimentally well verified for small and moderate Froude numbers $Fr\lesssim 0.3$ and frequencies much smaller than the resonance frequency, a regime that our theory should be able to describe. However, our solution (2.25) does not contain any linear scaling, except for the asymptotic limit of extremely small shaking frequencies $Fr\ll E$. In this limit, the surface elevation is determined solely by the balance between the centrifugal acceleration and the gravitational force. Figure 7(a) shows solution (2.25) (lines) together with the $\unicode[STIX]{x1D702}_{0}/R$ values (triangles) measured by Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) versus $Fr$ for different $E$ and $H$. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) used two glass cylinders of radii $R=5~\text{cm}$ and $6.5~\text{cm}$ partially filled with water of different depths. Weheliye et al.’s (Reference Weheliye, Yianneskis and Ducci2013) scaling law with $\unicode[STIX]{x1D6FC}_{0}=1.4$ is drawn in figure 7(a) as the dashed black line. Since all $E$ cases can be described by the potential model (except for the highest $Fr$, where the wave is becoming nonlinear), there is no linear dependency. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) have measured a series of resonance curves (compare with figure 4e) which appear to be linear in superposition. There is in fact only a linear dependency for the smallest $Fr<0.05$, which can be described by the flat disc solution $\unicode[STIX]{x1D702}_{0}/R=Fr$ (dash-dotted black line) contained in (2.32). For the sake of clarity, we have omitted the $H$ cases in figure 7(a); however, they fit to the potential model as well.

For Weheliye et al.’s (Reference Weheliye, Yianneskis and Ducci2013) set-up, solution (2.25) is not affected by damping rates of liquids with $\unicode[STIX]{x1D708}\lesssim 1000~\text{mm}^{2}~\text{s}^{-1}$ because all measurements were conducted in the non-resonant regime long before the first resonance is reached. As a result, equation (2.25) predicts the same surface elevations as the inviscid model of Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). This contradicts Ducci & Weheliye’s (Reference Ducci and Weheliye2014) assumption that the proportionality constant $\unicode[STIX]{x1D6FC}_{0}$ depends on the fluid viscosity. Based on $\unicode[STIX]{x1D702}_{0}/R$–$Fr$ measurements with oils of different viscosity, Ducci & Weheliye (Reference Ducci and Weheliye2014) found the following scaling law for $\unicode[STIX]{x1D6FC}_{0}$

(4.3)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{0}(\unicode[STIX]{x1D708})=\unicode[STIX]{x1D6FC}_{0w}\left({\displaystyle \frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D708}_{w}}}\right)^{-0.0256},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{0w}$ and $\unicode[STIX]{x1D708}_{w}$ denote the proportionality constant and the viscosity of water. The corresponding fitting curve is plotted in figure 1(b) of Ducci & Weheliye (Reference Ducci and Weheliye2014). In our opinion, this law does not properly describe the underlying physics for two reasons: first, it predicts a strong dependency on the viscosity for low-viscosity fluids close to water, while the dependency is very weak for high-viscosity liquids $\unicode[STIX]{x1D708}\gtrsim 100~\text{mm}^{2}~\text{s}^{-1}$. This is the exact opposite of what is predicted by our wave theory. Second, the decay of $\unicode[STIX]{x1D6FC}_{0}$ is extremely weak in the limit of high $\unicode[STIX]{x1D708}$ so that one obtains, e.g. for bitumen ($\unicode[STIX]{x1D708}\sim 1\times 10^{11}~\text{mm}^{2}~\text{s}^{-1}$) a constant of $\unicode[STIX]{x1D6FC}_{0}\approx 0.7$. Consequently, for bitumen, the law predicts sloshing with half the amplitude of water, an unlikely result.

Figure 7.

(a) Replot of the $E$ cases shown in (a) of Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013). Normalised surface elevations are plotted as a function of $Fr$. The markers show the measurements and the coloured lines represent the corresponding analytical predictions (2.25). In addition, the dashed black line marks the linear fit (4.2) with $\unicode[STIX]{x1D6FC}_{0}=1.4$ and the dash-dotted black line shows the disc solution $\unicode[STIX]{x1D702}_{0}/R=Fr$. (b) Predictions of critical $Fr$ numbers initiating the nonlinear out-of-phase flow transition as a function of $E$ for $H=1$ due to our wave theory-based prediction (4.6) (blue line) as well as the empirical law (4.5) presented in Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) (black line).

In our opinion, the concrete value of $\unicode[STIX]{x1D6FC}_{0}$ is mainly determined by the cutoff points chosen for the single measured $E$ and $H$ cases. If the highest $Fr$ measurement is removed for every $E$ case in figure 7(a), one would obtain a considerably lower fitting constant $\unicode[STIX]{x1D6FC}_{0}$. This can be seen from Ducci & Weheliye’s (Reference Ducci and Weheliye2014) figure 1(a) as well, where only small-$Fr$ cases are included, resulting in a reduced constant of $\unicode[STIX]{x1D6FC}_{0}=1.23$. Ducci & Weheliye (Reference Ducci and Weheliye2014) attribute this only to the higher viscosity of the employed liquid. However, except for the amplitude–viscosity dependency that we deem to be unphysical, equation (4.3) is a valuable estimate for the design of OSBs, but one has to know its limits.

Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) explain the deviation of the highest $Fr$ measurements from the fit with differences in surface waviness. The surface is flat for low $Fr$, but gets wavier with increasing $Fr$. This necessitates the use of an average of the local wave slopes, thereby introducing errors. Our interpretation of this observation is different. That the data points start to deviate from the fit as soon as the surface becomes wavy is in line with our theory of the disc–wave transition. In § 2.4 we have derived an expression for the critical $Fr$ number necessary to induce a certain waviness $\unicode[STIX]{x1D6F6}$

(4.4)$$\begin{eqnarray}Fr_{wave}(\unicode[STIX]{x1D6F6})={\displaystyle \frac{E\unicode[STIX]{x1D716}_{11}\tanh (\unicode[STIX]{x1D716}_{11}H)}{2(\unicode[STIX]{x1D716}_{11}^{2}-1)^{-1}\unicode[STIX]{x1D6F6}^{-1}+1}}.\end{eqnarray}$$

We found that all presented $E$ cases escape as soon as a waviness of $\unicode[STIX]{x1D6F6}=40\,\%$ is reached. All $E$ cases were conducted with $H=1$. Equation (4.4) predicts the critical $Fr_{wave}(40\,\%)=\{0.08,0.12,0.17,0.2,0.25,0.28\}$ for the plotted cases $E=\{0.14,0.21,0.3,0.36,0.44,0.5\}$. Figure 7(a) reveals that these $Fr_{wave}$ values coincide with the intersection points of the linear fit with the resonance curves. The evolution of a wavy-shaped surface is indeed a good indicator for a slowly emerging resonance here causing the nonlinear increases of the wave amplitudes. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) concluded correctly that $\unicode[STIX]{x1D6FC}_{0}$ is always greater than unity because it takes into account an extra inertial force, besides the centrifugal acceleration introduced by the shaker. This extra force is caused by the movement of the liquid itself; however, it is not constant and depends on the waviness and therefore also on $Fr$. Exactly this point was overlooked in the studies before.

To conclude, equation (4.4) can be used to determine the limits of the scaling law (4.2). Equation (4.2) provides a good estimate for surface elevations as long as $Fr\lesssim Fr_{wave}(\unicode[STIX]{x1D6F6})$. An empirical parameter in form of the waviness $\unicode[STIX]{x1D6F6}$ remains that depends on the fit. The smaller we choose the value of $\unicode[STIX]{x1D6F6}$ the more accurate (4.2) becomes. But physically there is no linear dependency in the $Fr$ regime considered and we recommend using (2.25) directly for more precise calculations.

Based on the preceding paragraphs, we can draw some conclusions regarding the out-of-phase flow frequently observed in OSBs. Such a flow is different from the linear phase shifts considered, e.g. by Bouvard et al. (Reference Bouvard, Herreman and Moisy2017) and Alpresa et al. (Reference Alpresa, Sherwin, Weinberg and van Reeuwijk2018a). The linear phase transition always arises at resonance and is smoothed out by internal damping forces. The phase transition observed in many OSB studies (Büchs et al. Reference Büchs, Maier, Milbradt and Zoels2000; Weheliye et al. Reference Weheliye, Yianneskis and Ducci2013; Ducci & Weheliye Reference Ducci and Weheliye2014; Klöckner et al. Reference Klöckner, Diederichs and Büchs2014; Thomas et al. Reference Thomas, Chakraborty, Berson, Shakeri and Sharp2017; Rodriguez et al. Reference Rodriguez, Micheletti and Ducci2018; Weheliye et al. Reference Weheliye, Cagney, Rodriguez, Micheletti and Ducci2018) is physically different and caused by nonlinear secondary flows. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) demonstrated the existence of two counter-rotating toroidal vortices close to the free surface, which are not predicted by our linear theory. With increasing $Fr$, these vortices move towards the walls and shrink until they disappear. This initiates the flow transition to the out-of-phase regime. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) observed that, simultaneously, the initially flat free surface becomes wavy, exactly the transformation which we have described using (4.4). Under the proposition that the out-of-phase transition always arises as soon as the flat disc is destabilised, we can also predict this nonlinear flow transition since the flat disc is always a linear solution captured by our theory. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) found the following empirical law for the critical $Fr$ number, above which the out-of-phase flow is expected to arise:

(4.5)$$\begin{eqnarray}\displaystyle Fr_{c}=\left\{\begin{array}{@{}ll@{}}{\displaystyle \frac{1}{2\unicode[STIX]{x1D6FC}_{0}}}H\sqrt{E}, & \text{if }H<2\sqrt{E},\\[7.0pt] {\displaystyle \frac{1}{\unicode[STIX]{x1D6FC}_{0}}}E, & \text{if }H>2\sqrt{E}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

On the basis of (4.4), we propose

(4.6)$$\begin{eqnarray}Fr_{c}={\displaystyle \frac{E\unicode[STIX]{x1D716}_{11}\tanh (\unicode[STIX]{x1D716}_{11}H)}{2(\unicode[STIX]{x1D716}_{11}^{2}-1)^{-1}\unicode[STIX]{x1D6F6}^{-1}+1}}\end{eqnarray}$$

as an alternative prediction. In direct comparison, it becomes apparent that (4.6) contains the predicted $H$ scaling of both cases in (4.5) as asymptotic limits: a linear scaling $Fr_{c}\sim H$ for $H\ll 1$ and a saturation $Fr_{c}=\text{const.}$ for $H\gg 1$. Hence, the $H$ dependence observed by Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) can be, at least partially, explained by the transition from the shallow water to the deep water wave regime. Both laws are graphically compared in figure 7(b) as a function of $E$ with a constant aspect ratio of $H=1$. Our criterion (4.6) also involves the empirical waviness $\unicode[STIX]{x1D6F6}$ (here set to $\unicode[STIX]{x1D6F6}=40\,\%$, comparable with $\unicode[STIX]{x1D6FC}_{0}=1.4$). It can be seen that (4.6) is very close to the established empirical law (4.5). Since the out-of-phase transition is not sharp and proceeds gradually, it can be well explained by the disc-to-wave transition alone. However, the ${\sim}\sqrt{E}$ scaling for small $H<2\sqrt{E}$ in (4.5) is not captured by (4.6). For larger $E$, equation (4.6) starts to overestimate critical $Fr$ numbers in comparison with (4.5). In this case, the flow transition is not caused by the motion of the free surface alone, which drives the toroidal vortices. Weheliye et al. (Reference Weheliye, Yianneskis and Ducci2013) found that for $2\sqrt{E}>H$ the vortices are steadily pushed towards the bottom wall. This initiates the out-of-phase transition in a way similar to the interaction with the sidewalls for $2\sqrt{E}<H$. In this case, waviness alone is not a sufficient criterion anymore and we cannot draw any further conclusions from our model. However, we can physically explain and support all other empirical scalings in (4.5) by relating the linear fitting constant $\unicode[STIX]{x1D6FC}_{0}$ to the specific waviness of the free surface.