Super-harmonically resonant swirling waves in longitudinally forced circular cylinders

Abstract

Resonant sloshing in circular cylinders was studied by Faltinsen et al. (J. Fluid Mech., vol. 804, 2016, pp. 608–645), whose theory was used to describe steady-state resonant waves due to a time-harmonic container's elliptic orbits. In the limit of longitudinal container motions, a symmetry breaking of the planar wave solution occurs, with clockwise and anti-clockwise swirling equally likely. In addition to this primary harmonic dynamics, previous experiments have unveiled that diverse super-harmonic dynamics are observable far from primary resonances. Among these, the so-called double-crest (DC) dynamics, first observed by Reclari et al. (Phys. Fluids, vol. 26, no. 5, 2014, 052104) for circular sloshing, is particularly relevant, as its manifestation is the most favoured by the spatial structure of the external driving. Following Bongarzone et al. (J. Fluid Mech., vol. 943, 2022, A28), in this work we develop a weakly nonlinear analysis to describe the system response to super-harmonic longitudinal forcing. The resulting system of amplitude equations predicts that a planar wave symmetry breaking via stable swirling may also occur under super-harmonic excitation. This finding is confirmed by our experimental observations, which identify three possible super-harmonic regimes, i.e. (i) stable planar DC waves, (ii) irregular motion and (iii) stable swirling DC waves, whose corresponding stability boundaries in the forcing frequency-amplitude plane quantitatively match the present theoretical estimates.

1. Introduction

Liquid sloshing related problems remain of great concern nowadays to many engineering fields. Depending on the type of external forcing and container shape, the free liquid surface can experience different types of dynamics, whose nature has a major importance in the design of, e.g. airplanes, rockets, spacecraft as well as road and ship tankers, since the sloshing motion may have a strong influence on their dynamic stability (Ibrahim 2005; Faltinsen & Timokha 2009). The case of resonant sloshing in upright circular cylinders represents one of the archetypal sloshing systems and it has indeed been extensively studied theoretically, experimentally and numerically.

In their work, Faltinsen, Lukovsky & Timokha (2016) thoroughly examine harmonically resonant sloshing dynamics in upright annular (circular) reservoirs. By applying the Narimanov–Moiseev multimodal sloshing theory (Narimanov 1957; Moiseev 1958; Dodge, Kana & Abramson 1965; Faltinsen 1974; Narimanov, Dokuchaev & Lukovsky 1977; Lukovsky 1990; Lukovsky & Timokha 2011, 2015; Takahara & Kimura 2012; Lukovsky 2015), capable of accurately describing the nonlinear wave dynamics near primary harmonic resonances and in the absence of secondary resonances (Faltinsen, Rognebakke & Timokha 2005; Faltinsen et al. 2016; Raynovskyy & Timokha 2018a, 2020), i.e. for a non-dimensional fluid depth $H\gtrsim 1.05$, they derived the response curves for planar elliptic-type tank excitation. In the two limit cases, system responses to longitudinal and circular tank motions were retrieved.

Circular sloshing is widely used in biological and chemical industrial applications such as small- and large-scale bioreactors for bacterial and cellular cultures (McDaniel & Bailey 1969; Wurm 2004), where the liquid motion prevents the sedimentation of suspended cells in the liquid medium and allows for a homogenized concentration of dissolved oxygen and nutrients. For these reasons, a strong interest in the gas exchange and mixing processes taking place in these devices has emerged over the last decades (Büchs et al. 2000a,b; Büchs 2001; Maier, Losen & Büchs 2004; Muller et al. 2005; Micheletti et al. 2006; Zhang et al. 2009; Tissot et al. 2010; Tan, Eberhard & Büchs 2011; Tissot et al. 2011; Klöckner & Büchs 2012).

Reclari (2013) and Reclari et al. (2014), among others (see also Hutton 1964; Bouvard, Herreman & Moisy 2017; Moisy, Bouvard & Herreman 2018; Horstmann, Herreman & Weier 2020; Horstmann et al. 2021), experimentally characterized in great detail the hydrodynamics of orbitally shaken circular cylinders, which represent the typical shape of lab-scale bioreactors. In addition to the primary harmonic system response via single-crest (SC) wave dynamics, different multiple-crest wave patterns were observed. Among these, the super-harmonic double-crest (DC) wave dynamics, as labelled by Reclari et al. (2014), is particularly relevant, as it appears to be the most stable and the one that displays the largest nonlinear amplitude response, that may eventually lead to wave breaking occurring far from harmonic resonances and even at moderately low forcing amplitudes. Its manifestation is indeed naturally favoured by the spatial structure, i.e. by the temporal and azimuthal periodicities, of the external driving and, therefore, its understanding and prediction can be important for practical application as in the design of bioreactors.

The analysis outlined in Bongarzone, Guido & Gallaire (2022a) was precisely dedicated to the development of an inviscid weakly nonlinear (WNL) analysis, which was seen to successfully capture nonlinear effects for this subtle additive and multiplicative resonance governing the super-harmonic DC swirling and which well matched the experimental findings of Reclari (2013) and Reclari et al. (2014).

Nonetheless, the applicability of the aforementioned analysis is limited to circular sloshing, whereas the emergence of super-harmonic DC dynamics is in principle expected for any elliptic-type container excitation and, therefore, for longitudinal forcing as well.

The latter forcing condition has been analytically and experimentally studied for decades (Hutton 1963; Abramson 1966; Chu 1968) and it is of interest from the perspective of hydrodynamic instabilities due to the occurrence of hysteretic symmetry-breaking conditions (Miles 1984a,b). With regards to circular cylindrical containers, particularly relevant are the experimental studies by Abramson, Chu & Kana (1966), Royon-Lebeaud, Hopfinger & Cartellier (2007) and Hopfinger & Baumbach (2009), who detected the stability bounds between harmonic planar, swirling and irregular waves and whose estimates were later used by Faltinsen et al. (2016) to validate their theoretical analysis. However, these works were mostly focused on the investigation of system responses in the neighbourhood of harmonic resonances, whereas, with the exception of Reclari et al. (2014) and Bongarzone et al. (2022a) in the context of circular sloshing, the literature seems to lack comprehensive experimental and theoretical studies dealing with the most relevant secondary super-harmonic resonances (by super harmonic, we mean here a wave of a certain frequency $\omega$ emerging from an excitation at $\varOmega = \omega /2$, with $\varOmega$ the driving angular frequency), i.e. far from primary ones, under longitudinal or, more generally, elliptical container excitation.

In this work we take a first step in this direction by applying to longitudinal planar forcing the analysis formalized by Bongarzone et al. (2022a) for circular container motions. In the spirit of the multiple time-scale method, we develop a WNL model leading to a system of two amplitude equations, which, via thorough comparison with dedicated lab-scale experiments, is proven capable of describing satisfactorily the steady-state system response to super-harmonic longitudinal forcing and, particularly, of detecting the various possible dynamical regimes.

The manuscript is organized as follows. The flow configuration and governing equations are given in § 2. In § 3 we briefly introduce the classical linear potential model and briefly describe the numerical method employed in this work. By analogy with Bongarzone et al. (2022a), in § 4 we first tackle the simpler case of the harmonic SC wave. The WNL system of amplitude equations governing the DC wave dynamics under super-harmonic longitudinal forcing, which represents the core of this study, is then formalized in § 5. The experimental apparatus, procedure and findings are described in § 6, where a thorough quantitative comparison with the present theoretical estimates is carried out. Final comments and conclusions are outlined in § 7. Lastly, Appendix B complements the theoretical model by briefly showing how a straightforward extension of the present analysis to a generic container's elliptic orbits can be readily obtained without any further calculation, hence paving the way for further analyses and experimental investigations.

2. Flow configuration and governing equations

We consider a cylindrical container of diameter $D=2R$ filled to a depth $h$ with a liquid of density $\rho$. The air–liquid surface tension is denoted by $\gamma$, whereas the gravity acceleration is denoted by $g$. Here $O'\boldsymbol {e}'_x\boldsymbol {e}'_y\boldsymbol {e}'_z$ is the Cartesian inertial reference frame, while $O\boldsymbol {e}_x\boldsymbol {e}_y\boldsymbol {e}_z$ is the Cartesian reference frame moving with the container. The origin of the moving cylindrical reference frame $(r,\theta,z)$ is placed at the container revolution axis and, specifically, at the unperturbed liquid height, $z=0$ (see figure 1). A longitudinal shaking in the horizontal plane, e.g. along the $x$ axis, can be represented by the following equations describing the motion velocity of the container axis intersection with the $z=0$ plane, parametrized in polar coordinates ($r$, $\theta$):

(2.1)\begin{equation} \dot{\boldsymbol{X}}_0 = \begin{cases} - \bar{a}_x{\bar{\varOmega}}\sin{({\bar{\varOmega} }t)}\cos\theta \boldsymbol{e}_r,\\ \bar{a}_x{\bar{\varOmega}}\sin{({\bar{\varOmega}} t)}\sin\theta \boldsymbol{e}_{\theta}. \end{cases} \end{equation}

Here $\bar {a}_x$ is the dimensional forcing amplitude and $\bar {\varOmega }$ is the dimensional driving angular frequency. In the potential flow limit, the liquid motion within the moving container is governed by the Laplace equation, subjected to the homogeneous no-penetration condition at the solid lateral wall and bottom,

(2.2a,b)\begin{equation} {\rm \Delta}\varPhi=0,\quad \boldsymbol{\nabla}\varPhi\boldsymbol{\cdot}\boldsymbol{n}= \boldsymbol{0}, \end{equation}

and by the dynamic and kinematic boundary conditions at the free surface $z=\eta (r,\theta )$ (Ibrahim 2005; Faltinsen & Timokha 2009),

(2.3a)\begin{gather} \frac{\partial\varPhi}{\partial t}+\frac{1}{2} \boldsymbol{\nabla}\varPhi\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi+\eta- \frac{\kappa(\eta)}{Bo}=rf\cos{(\varOmega t)} \cos{\theta}, \end{gather}

(2.3b)\begin{gather}\frac{\partial \eta}{\partial t}+\boldsymbol{\nabla}\varPhi\boldsymbol{\cdot}\boldsymbol{\nabla}\eta-\frac{\partial\varPhi}{\partial z}=0, \end{gather}

which have been made non-dimensional by using the container's characteristic length $R$, the velocity $\sqrt {gR}$ and the time scale $\sqrt {R/g}$. In (2.3a), $\kappa (\eta )$ denotes the fully nonlinear curvature, while $Bo=\rho gR^2/\gamma$ is the Bond number. As soon as the Bond number is sufficiently large, i.e. $Bo\sim 10^3$ (Bouvard et al. 2017), surface tension effects are almost negligible (fully negligible for $Bo \gtrsim 10^4$, except in the neighbourhood of the contact line Faltinsen et al. 2016). In the following, we assume large Bond numbers and, accordingly, the curvature term in (2.3a) is neglected. The non-dimensional driving acceleration along the $x$ axis reads $f=a_x\varOmega ^2$, with $a_x=\bar {a}_x/R$ and $\varOmega ={\bar {\varOmega }}/\sqrt {g/R}$. Lastly, the non-dimensional fluid depth is $H=h/R$.

Figure 1. Sketch of a cylindrical container of diameter $D=2R$ and filled to a depth $h$. The gravity acceleration is denoted by $g$. Here $O'\boldsymbol {e}'_x\boldsymbol {e}'_y\boldsymbol {e}'_z$ is the Cartesian inertial reference frame, while $O\boldsymbol {e}_x\boldsymbol {e}_y\boldsymbol {e}_z$ is the Cartesian reference frame moving with the container. The origin of the moving cylindrical reference frame $(r,\theta,z)$ is placed at the container revolution axis and, specifically, at the unperturbed liquid height, $z=0$. The perturbed free-surface and contact line elevation are denoted by $\eta$ and $\delta$, respectively. Here $\bar {a}_x$ is the amplitudes of the longitudinal periodic forcing of the driving angular frequency $\bar {\varOmega }$.

3. Linear potential model

Far from resonances and in the limit of small forcing amplitudes, the linear theory is expected to provide a good approximation of the harmonic system response. Let us consider small perturbations of a base state $\boldsymbol {q}_0$,

(3.1)\begin{equation} \boldsymbol{q}(r,\theta,z,t)= \{\varPhi(r,\theta,z,t),\eta(r,\theta,t)\}^T= \boldsymbol{q}_0+\epsilon \boldsymbol{q}'= \epsilon\{\varPhi',\eta'\}^T+{O} (\epsilon^2), \end{equation}

together with the assumption of small driving forcing amplitudes of order ${O}(\epsilon )$, i.e. $f=\epsilon F$, with $\epsilon$ a small parameter $\epsilon \ll 1$ and with the auxiliary variable $F$ of order ${O}(1)$.

In the following, we assume that the dynamic contact line freely slides along the lateral wall with a constant slope while keeping a contact angle equal to a static value of $90^{\circ }$. The latter hypothesis implicitly assumes the absence of a static meniscus so that the base-state configuration is $\boldsymbol {q}_0=\{\varPhi _0,\eta _0\}^T=\boldsymbol {0}$, i.e. the fluid is at rest with a flat static interface.

At order $\epsilon$, (2.2a,b) and (2.3b) reduce to a forced linear system, whose matrix compact form reads

(3.2)\begin{equation} (\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}) \boldsymbol{q}'=\boldsymbol{\mathcal{F}}', \end{equation}

with $\boldsymbol {\mathcal {F}}'=F\hat {\boldsymbol {\mathcal {F}}} (\frac {1}{2}\, {\rm e}^{\text {i}(\varOmega t-\theta )}+\frac {1}{2}\, {\rm e}^{\text {i}(\varOmega t+\theta )})+{\rm c.c.}$, $\hat {\boldsymbol {\mathcal {F}}}=\{0,{r}/{2}\}^T$ and

(3.3a,b)\begin{equation} \boldsymbol{\mathsf{B}} = \begin{pmatrix} 0 & 0 \\ I_{\eta} & 0\\ \end{pmatrix},\quad \boldsymbol{\mathsf{A}} = \begin{pmatrix} {\rm \Delta} & 0 \\ 0 & -I_{\eta}\\ \end{pmatrix}, \end{equation}

where c.c. stands for complex conjugate and $I_{\eta }$ is the identity matrix associated with the interface $\eta$. Note that the kinematic condition does not explicitly appear in (3.3a,b), but it is enforced as a boundary condition at the interface (Viola, Brun & Gallaire 2018). We then seek a standing wave solution in the form

(3.4)\begin{equation} \boldsymbol{q}'(r,\theta,z,t)=F\hat{\boldsymbol{q}} (r,z)(\tfrac{1}{2}\, {\rm e}^{\text{i} (\varOmega t-\theta)}+\tfrac{1}{2}\, {\rm e}^{\text{i} (\varOmega t+\theta)})+{\rm c.c.}, \end{equation}

where $\hat {\boldsymbol {q}}$ is straightforwardly computed by solving the system

(3.5)\begin{equation} (\text{i}\varOmega\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m=1}}) \hat{\boldsymbol{q}}=\hat{\boldsymbol{\mathcal{F}}}. \end{equation}

Note that, due to the normal mode ansatz (3.4), the linear operator $\boldsymbol{\mathsf{A}}_m$ depends on the azimuthal wavenumber $m$, here $m=1$. Even though an exact analytical solution to (3.5) can be readily obtained via a Bessel–Fourier series representation, in this work, as in Bongarzone et al. (2022a), we opt for a numerical scheme based on a discretization technique, where linear operators $\boldsymbol{\mathsf{B}}$ and $\boldsymbol{\mathsf{A}}_{m}$ are discretized in space by means of a Chebyshev pseudo-spectral collocation method with a two-dimensional mapping implemented in Matlab, which is analogous to that described by Viola et al. (2018) and Bongarzone, Viola & Gallaire (2021b). The numerical scheme requires explicit boundary conditions at $r=0$ in order to regularize the problem on the revolution axis ($r=0$), i.e.

(3.6a)\begin{gather} m=0:\ \frac{\partial\hat{\eta}}{\partial r}=\frac{\partial\hat{\varPhi}}{\partial r}=0, \end{gather}

(3.6b)\begin{gather}{m}\ge1:\ \hat{\eta}=\hat{\varPhi}=0. \end{gather}

The numerical convergence of the results presented throughout the paper is achieved using a computational grid $N_r = N_z \le 40$, with $N_r$ and $N_z$ the number of radial and axial collocation grid points, respectively. Due to the low computational cost, we used $N_r = N_z = 60$.

We recall the well-known dispersion relation for inviscid gravity waves (Lamb 1993),

(3.7)\begin{equation} \omega_{mn}^2=k_{mn}\tanh{(k_{mn}H)}, \end{equation}

where the wavenumber $k_{mn}$ is given by the nth root of the first derivative of the mth-order Bessel function of the first kind satisfying $J'_{m}(k_{mn})=0$. By denoting the eigenvector associated with the natural frequency $\omega _{mn}$ as $\hat {\boldsymbol {q}}_{mn}$, the solution of the homogeneous version of (3.5) for $\varOmega =\omega _{mn}$, it is useful for the rest of the analysis to note that owing to the symmetries of the problem, the system admits the following invariant transformation:

(3.8)\begin{equation} (\hat{\boldsymbol{q}}_{mn},+m,\text{i}\omega_{mn}) \longrightarrow(\hat{\boldsymbol{q}}_{mn},-m,\text{i}\omega_{mn}). \end{equation}

Such an invariance suggests that the spatial structure, $\hat {\boldsymbol {q}}(r,z)$, of the system response to an external forcing with temporal and azimuthal periodicity $(\varOmega,m)$ is the same as that computed for $(\varOmega,-m)$, so that the linear solution form (3.4) holds true.

4. Harmonic SC resonance

With the aim to derive a WNL system of amplitude equations governing the super-harmonic DC dynamics under longitudinal excitation, we first tackle the simpler problem of harmonic SC waves. We look for a third-order asymptotic solution of the system

(4.1)\begin{equation} \boldsymbol{q}=\{\varPhi,\eta\}^T=\epsilon \boldsymbol{q}_1+\epsilon^2\boldsymbol{q}_2+ \epsilon^3\boldsymbol{q}_3+{O}(\epsilon^4), \end{equation}

where the zero-order solution, $\boldsymbol {q}_0=\boldsymbol {0}$, associated with the rest state, is omitted.

With regards to SC waves and, specifically, to the harmonic response at a driving frequency close to the natural frequency of one of the non-axisymmetric $m=\pm 1$ modes, $\omega _{1n}$, we assume here a small forcing amplitude of order $\epsilon ^3$. This assumption is justified by the fact that close to resonance, $\varOmega \approx \omega _{1n}$, and in the absence of dissipation, even a small forcing will induce a large system response. Hence, the analysis is expected to hold for $\varOmega =\omega _{1n}+\lambda$, where $\lambda$ is a small detuning parameter assumed of order $\epsilon ^2$. In the spirit of the multiple-scale technique, we introduce the slow time scale $T_2 = \epsilon ^2t$, with $t$ being the fast time scale. Hence, the scalings

(4.2a–c)\begin{equation} f=\epsilon^3 F,\quad \varOmega=\omega_{1n}+\epsilon^2\varLambda, \quad T_2=\epsilon^2 t \end{equation}

are assumed, with the auxiliary parameters, $F$ and $\varLambda$, of order ${O}(1)$.

Given the azimuthal periodicity of the external forcing, i.e. $m=\pm 1$, we assume a leading-order solution as the sum of two counter-propagating travelling waves,

(4.3)\begin{equation} \boldsymbol{q}_1=A_1(T_2)\hat{\boldsymbol{q}}_{1}^{A_1}\, {\rm e}^{\text{i}(\omega_{1n}t-\theta)}+B_1(T_2) \hat{\boldsymbol{q}}_{1}^{B_1}\, {\rm e}^{\text{i} (\omega_{1n}t+\theta)}+{\rm c.c.}, \end{equation}

where $\hat {\boldsymbol {q}}_1^{A_1}=\hat {\boldsymbol {q}}_1^{B_1}$ (owing to (3.8)) is the eigenmode computed by solving (3.5) for its homogeneous solution at $\varOmega =\omega _{1n}$, where $\omega _{1n}$ is given by (3.7). The complex amplitudes $A_1$ and $B_1$, functions of the slow time scale $T_2$ and still undetermined at this stage of the expansion, describe the slow time amplitude modulation of the two oscillating waves and must be determined at a higher order.

By pursuing the expansion to the second order, one obtains a linear system forced by combinations of the first-order solutions. These forcing terms are proportional to $A_1^2$ and $B_1^2$ (second harmonics), to $|A_1|^2$ and $|B_1|^2$ (steady and axisymmetric mean flow corrections) and to $A_1B_1$ and $A_1\bar {B}_1$ (cross-quadratic interactions),

(4.4)\begin{align} &(\partial_t\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m}}) \boldsymbol{q}_2=\mathcal{F}_2=(|A_1|^2 \hat{\boldsymbol{\mathcal{F}}}_2^{A_1\bar{A}_1}+|B_1|^2_1 \hat{\boldsymbol{\mathcal{F}}}_2^{B_1\bar{B}_1})\nonumber\\ &\quad +(A_1^2\hat{\boldsymbol{\mathcal{F}}}_{2}^{A_1A_1}\, {\rm e}^{\text{i}2(\omega_{1n}t-\theta)}+B_1^2 \hat{\boldsymbol{\mathcal{F}}}_{2}^{B_1B_1}\, {\rm e}^{\text{i}2 (\omega_{1n}t+\theta)}+{\rm c.c.})\nonumber\\ &\quad +(A_1B_1\hat{\boldsymbol{\mathcal{F}}}_{2}^{A_1B_1}\, {\rm e}^{\text{i}2\omega_{1n}t}+A_1\bar{B}_1 \hat{\boldsymbol{\mathcal{F}}}_{2}^{A_1\bar{B}_1}\, {\rm e}^{-\text{i}2\theta}+{\rm c.c.}). \end{align}

Thus, we seek for a second-order solution of the form

(4.5)\begin{align} \boldsymbol{q}_2&=|A_1|^2\hat{\boldsymbol{q}}_2^{A_1\bar{A}_1} +|B_1|^2\hat{\boldsymbol{q}}_2^{B_1\bar{B}_1}+ (A_1^2\hat{\boldsymbol{q}}_{2}^{A_1A_1}\, {\rm e}^{\text{i}2 (\omega_{1n}t-\theta)}+B_1^2\hat{\boldsymbol{q}}_{2}^{B_1B_1}\, {\rm e}^{\text{i}2(\omega_{1n}t+\theta)}+{\rm c.c.})\nonumber\\ &\quad +(A_1B_1\hat{\boldsymbol{q}}_{2}^{A_1B_1}\, {\rm e}^{\text{i}2\omega_{1n}t}+A_1\bar{B}_1 \hat{\boldsymbol{q}}_{2}^{A_1\bar{B}_1}\, {\rm e}^{-\text{i} 2\theta}+{\rm c.c.}). \end{align}

Given the invariant transformation (3.8), only some of these second-order responses need to be computed explicitly, as, e.g. $\hat {\boldsymbol {q}}_2^{A_1\bar {A}_1}= \hat {\boldsymbol {q}}_2^{B_1\bar {B}_1}$ and $\hat {\boldsymbol {q}}_2^{A_1A_1}=\hat {\boldsymbol {q}}_2^{B_1B_1}$.

We now move forward to the $\epsilon ^3$-order problem, which is once again a linear problem forced by combinations of the first-order (4.3) and second-order solutions (4.5), produced by third-order nonlinearities such as $(\boldsymbol{\nabla} \varPhi _1\boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi _2+\boldsymbol{\nabla} \varPhi _2 \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi _1)/2$ in the dynamic condition or $\boldsymbol{\nabla} \varPhi _1\boldsymbol{\cdot} \boldsymbol{\nabla} \eta _2+\boldsymbol{\nabla} \varPhi _2\boldsymbol{\cdot} \boldsymbol{\nabla} \eta _1$ in the kinematic equation, as well as by the slow time $T_2$ derivative of the leading-order solution and by the external forcing, which was assumed of order $\epsilon ^3$,

(4.6)\begin{align} &(\partial_t\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m}}) \boldsymbol{q}_3=\boldsymbol{\mathcal{F}}_3={-}\frac{\partial A_1}{\partial T_2} \boldsymbol{\mathsf{B}} \hat{\boldsymbol{q}}_1^{A_1}\, {\rm e}^{\text{i} (\omega_{1n}t-\theta)}-\frac{\partial B_1}{\partial T_2} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{B_1}\, {\rm e}^{\text{i}(\omega_{1n}t+\theta)}\nonumber\\ &\quad + |A_1|^2A_1\hat{\boldsymbol{\mathcal{F}}}_3^{|A_1|^2A_1}\, {\rm e}^{\text{i}(\omega_{1n}t-\theta)}+ |B_1|^2B_1 \hat{\boldsymbol{\mathcal{F}}}_3^{|B_1|^2B_1}\, {\rm e}^{\text{i}(\omega_{1n}t+\theta)}\nonumber\\ &\quad + |B_1|^2A_1 \hat{\boldsymbol{\mathcal{F}}}_3^{|B_1|^2A_1}\, {\rm e}^{\text{i}(\omega_{1n}t-\theta)}+ |A_1|^2B_1 \hat{\boldsymbol{\mathcal{F}}}_3^{|A_1|^2B_1}\, {\rm e}^{\text{i}(\omega_{1n}t+\theta)}\nonumber\\ &\quad+\frac{1}{2}F\hat{\boldsymbol{\mathcal{F}}}_3^F {\rm e}^{\text{i}(\omega_{1n}t-\theta)}\, {\rm e}^{\text{i}\varLambda T_2}+\frac{1}{2}F\hat{\boldsymbol{\mathcal{F}}}_3^F {\rm e}^{\text{i}(\omega_{1n}t+\theta)}\, {\rm e}^{\text{i}\varLambda T_2} \nonumber\\ &\quad +\text{N.R.T.}+{\rm c.c.}, \end{align}

with $\hat {\boldsymbol {\mathcal {F}}}_3^F=\{0,r/2\}^T$ and where $\text {N.R.T.}$ stands for non-resonating terms. These terms are not strictly relevant for further analysis and can therefore be neglected. Amplitude equations for $A_1$ and $B_1$ are obtained by requiring that secular terms do not appear in the solution to (4.6), where secularity results from all resonant forcing terms in $\boldsymbol {\mathcal {F}}_3$ (see Appendix D of Bongarzone et al. (2022a) for its explicit expression), i.e. all terms sharing the same frequency and wavenumber of $\boldsymbol {q}_1$, e.g. $(\omega,m)=(\omega _{1n},\pm 1)$, and in effect all terms explicitly written in (4.6). It follows that a compatibility condition must be enforced through the Fredholm alternative (Friedrichs 2012; Olver 2014), which imposes the amplitudes $A=\epsilon A_1e^{-\text {i}\lambda t}$ and $B=\epsilon B_1e^{-\text {i}\lambda t}$ to obey the following normal form:

(4.7a)\begin{gather} \frac{{\rm d} A}{{\rm d} t}={-}\text{i}\lambda A + \text{i}\frac{\mu_{{SC}}}{2}f + \text{i}\nu_{{SC}} |A|^2A +\text{i}\xi_{{SC}}|B|^2A, \end{gather}

(4.7b)\begin{gather}\frac{{\rm d} B}{{\rm d} t}={-}\text{i}\lambda B + \text{i}\frac{\mu_{{SC}}}{2}f + \text{i}\nu_{{SC}} |B|^2B + \text{i}\xi_{{SC}}|A|^2B. \end{gather}

Here the physical time $t=T_2/\epsilon ^2$ has been reintroduced and where the forcing amplitude and detuning parameter are recast in terms of their corresponding physical values, $f=\epsilon ^3 F$ and $\lambda =\epsilon ^2\varLambda =\varOmega -\omega _{1n}$, so as to eliminate the small implicit parameter $\epsilon$ (Bongarzone et al. 2021a, 2022b). The subscript $_{SC}$ stands for single crest. The various normal form coefficients, which turn out to be real-valued quantities due to the absence of dissipation, are computed as scalar products between the adjoint mode, $\hat {\boldsymbol {q}}_1^{A_1 {\dagger} }=\hat {\boldsymbol {q}}_1^{B_1 {\dagger} }$, associated with $\hat {\boldsymbol {q}}_1^{A_1}=\hat {\boldsymbol {q}}_1^{B_1}$, and the third-order resonant forcing terms (see Appendix A and Bongarzone et al. (2022a) for further details).

Once stable stationary solutions are computed, $A$ and $B$ are replaced in expressions (4.3) and (4.5) and the total harmonic SC wave solution is reconstructed as

(4.8)\begin{equation} \boldsymbol{q}_{SC}=\{\varPhi,\eta\}^T= \epsilon\boldsymbol{q}_1+\epsilon^2\boldsymbol{q}_2. \end{equation}

To this end, it is first convenient to express (4.7a) and (4.7b) in polar coordinates, i.e. by defining $A=|A|{\rm e}^{\text {i}\varPhi _A}$ and $B=|B|{\rm e}^{\text {i}\varPhi _B}$, and then to introduce the following change of variables, $|a|=|A|+|B|$ and $|b|=|A|-|B|$. By looking for periodic solutions with stationary amplitudes $|A|,\,|B|\ne 0$, one can sum and subtract (4.7a) and (4.7b), hence obtaining

(4.9a)\begin{gather} f=a_x\varOmega^2={\pm} |a|\left(\lambda-\left(\frac{\nu_{{SC}}+ \xi_{SC}}{4}\right)|a|^2-\left(\frac{3\nu_{{SC}}- \xi_{SC}}{4}\right)|b|^2\right)\frac{1}{\mu_{{SC}}}, \end{gather}

(4.9b)\begin{gather}0=|b|\left(\lambda-\left(\frac{\nu_{{SC}}+\xi_{{SC}}}{4} \right)|b|^2-\left(\frac{3\nu_{{SC}}-\xi_{SC}}{4}\right)|a|^2\right). \end{gather}

As expected, (4.9b) suggests that two possible solutions exist. The planar (or standing) wave solution is retrieved for

(4.10a)\begin{gather} |b|=|A|-|B|=0\rightarrow |A|=|B|, \end{gather}

(4.10b)\begin{gather}a_x\varOmega^2={\pm} |a|\left(\lambda-\left(\frac{\nu_{{SC}}+ \xi_{{SC}}}{4}\right)|a|^2\right)\frac{1}{\mu_{{SC}}}, \end{gather}

whereas the swirling wave solution is found when $|b|\ne 0$ and

(4.11a)\begin{gather} |b|^2=\left(\lambda-\left(\frac{3\nu_{{SC}}-\xi_{{SC}}}{4} \right)|a|^2\right)\left(\frac{4}{\nu_{{SC}}+\xi_{{SC}}}\right), \end{gather}

(4.11b)\begin{gather}a_x\varOmega^2={\pm} 2|a|\left(\frac{\xi_{{SC}}-\nu_{{SC}}}{\nu_{{SC}}+ \xi_{{SC}}}\right)\left(\lambda-\nu_{{SC}}|a|^2\right) \frac{1}{\mu_{{SC}}}. \end{gather}

The various branches prescribed by (4.10b), (4.11a) and (4.11b) for $|a|$ and $|b|$ as a function of $\tau =\varOmega /\omega _{1n}$ and at a fixed non-dimensional shaking amplitude $a_x$ are here computed by means of the Matlab function fimplicit.

From (4.7a) and (4.7b) expressed in polar coordinates, one finds that the stationary module equations read $\mu _{{SC}}f\sin {\varPhi _A}/2=0$ and $\mu _{{SC}}f\sin {\varPhi _B}/2=0$, hence implying $\sin {\varPhi _A}=\sin {\varPhi _B}=0$. We therefore note that four possible combinations of stationary phases, $\varPhi _A$ and $\varPhi _B\in [0,2{\rm \pi} ]$, are in principle admitted, i.e. (i) $\varPhi _A=\varPhi _B=0$, (ii) $\varPhi _A=\varPhi _B={\rm \pi}$, (iii) $\varPhi _A=0$, $\varPhi _B={\rm \pi}$ and (iv) $\varPhi _A={\rm \pi}$, $\varPhi _B=0$. However, (iii) and (iv) are totally equivalent to (i) and (ii), respectively, with amplitudes $|a|\rightarrow |b|$ and $|b|\rightarrow |a|$. Therefore, only combinations (i) $\varPhi _A=\varPhi _B=\varPhi =0$ and (ii) $\varPhi ={\rm \pi}$, which produce the $\pm$ sign in (4.9a), are retained.

4.1. Comparison with existing experiments and theoretical predictions

In figure 2 we reproduce figure 8 of Faltinsen et al. (2016), which shows the estimates of bounds between the frequency ranges where harmonic planar, irregular and swirling waves occur. The outcomes of the present analysis are consistent with those of Faltinsen et al. (2016) and with the experimental measurements by Royon-Lebeaud et al. (2007). The values of the normal form coefficients $\mu _{{SC}}$, $\nu _{{SC}}$ and $\xi _{{SC}}$ reported in table 1 of Appendix A confirm that the stability boundaries vary weakly with the liquid depth, as stated by Faltinsen et al. (2016) for non-dimensional fluid depths $H\gtrsim 1.05$, but strongly depend on the forcing amplitude, with the frequency range for irregular and swirling waves widening for increasing forcing amplitudes. In this context, irregular means that both the planar and the swirling wave solutions are unstable, hence, one could expect irregular and chaotic patterns with a switching between planar and swirling motion. The green shaded region corresponds to stable SC swirling waves, while the light purple shaded region corresponds to the multi-solution regime, where both stable swirling SC and planar SC wave motions are possible depending on the initial transient, i.e. on the initial conditions, as typical of hysteretic systems.

Figure 2. Estimates of bounds, in the $(\varOmega /\omega _{11},a_x)$ plane, between the frequency ranges where planar, irregular and swirling waves occur when the container undergoes a longitudinal and harmonic motion. Filled markers: experiments by Royon-Lebeaud et al. (2007). Black dashed lines: theoretical prediction by Faltinsen et al. (2016), whose theoretical curves have been reproduced here by manually sampling those reported in their original figure 8(a).

Table 1. Value of the normal form coefficients appearing in (4.7a) and (4.7b) (SC) and in (5.9a) and (5.9b) (DC) computed at different fluid depths $H=h/R$ and associated with mode $(m,n)=(1,1)$. Note that in (5.9a) and (5.9b), $\chi _{{DC}}=\chi _{{-}}+\chi _{{+}}$.

In figures 3(a) and 3(b) the non-dimensional maximum steady-state wave elevation, computed by reconstructing the total flow solution in accordance with (4.8), is compared with the theoretical estimations by Faltinsen et al. (2016) (black dashed lines) from their figure 8 and with the corresponding experimental measurements by Royon-Lebeaud et al. (2007) (coloured filled markers). The agreement between the present model and experiments is fairly good and consistent with predictions by Faltinsen et al. (2016). The larger disagreement between theory and experiments at smaller forcing amplitudes was tentatively attributed by Faltinsen et al. (2016) to the fact that the actual elevation of these wave amplitudes was approximately $1\ \text {mm}$ and may therefore be more difficult to measure with sufficient accuracy.

Figure 3. Non-dimensional maximum steady-state wave elevation, $\max _{t,\theta =0,{\rm \pi} /2}\,\eta$ (the maximum is taken from values at two probes located at $(x,y)=(0.875,0)$ and $(0,0.875)$) versus the forcing frequency $\varOmega /\omega _{11}$ and for different $x$-longitudinal shaking amplitudes, $a_x$: (a) $0.0033$, $0.0066$, $0.0133$ and $0.0266$; (b) $0.023$ and $0.045$. Markers are associated with two experimental series by Royon-Lebeaud et al. (2007) (experimental data from their original figures 2 (now (a)) and 7 (now (b)). Filled circles correspond to measurements done for the planar regime, whereas filled squares indicate swirling. The black dashed lines represent the stable branches predicted by Faltinsen et al. (2016). Their curves have been carefully reproduced here by manually sampling those reported in their original figure 10 in the range of frequency available, i.e. $\varOmega /\omega _{11}\in [0.7,1.2]$. Coloured solid lines correspond to the present theoretical predictions for stable branches.

A comparable mismatch is here retrieved. As a side comment, we note that, within the present inviscid framework, the lower left stable planar branch, $\varOmega /\omega _{11}<1$, is obtained for a phase $\varPhi =0$, which implies a fluid motion in phase with the container motion, whereas the lower right planar branch, $\varOmega /\omega _{11}>1$, has a phase $\varPhi ={\rm \pi}$, hence implying a phase opposition. The stable swirling branch is characterized by $\varPhi =0$. This is consistent with previous studies (Royon-Lebeaud et al. 2007).

4.2. Discussion: present analysis vs the Narimanov–Moiseev multimodal theory

In this section the present model for harmonic resonances has been compared with the Narimanov–Moiseev multimodal theory employed by Faltinsen et al. (2016) and has been shown to provide very consistent and close predictions. Before moving to the next section, it is therefore worth pointing out the methodological analogies and differences as well as the pros and cons of the two approaches.

Adopting a variational formulation (Miles 1976) and assuming an incompressible and irrotational flow, the multimodal method reduces the hydrodynamic sloshing system to a modal system of nonlinearly coupled ordinary differential equations written in terms of the so-called generalized coordinates (Faltinsen & Timokha 2009). This projection step uses a Fourier-type representation of the time-dependent surface elevation and potential velocity field. Because the resulting coefficients in the equation system are derived analytically and only ordinary differential equations must be numerically time integrated, numerical errors are negligible with small CPU time relative to that for computational fluid dynamics methods based on the governing equations of the full hydrodynamic sloshing system.

For theoretical modelling purposes, postulating proper asymptotic relations between these generalized coordinates simplifies the system to a WNL form. Specifically in the case of harmonic resonances, the Narimanov–Moiseev asymptotic relations assume a leading dynamics of order $\sim \epsilon$, a frequency detuning $\sim \epsilon ^2$, a forcing amplitude $\sim \epsilon ^3$ and a slow time scale $\sim \epsilon ^2$. As surface tension effects are neglected and the contact line is assumed to freely sleep along the sidewall with a constant and zero slope, the Fourier basis (Bessel functions for circular cylinders) also coincides with the actual natural sloshing modes.

Under these assumptions, although we do not go through this initial projection step, the present model and the Narimanov–Moiseev multimodal theory are essentially equivalent. The $\epsilon$-order leading dynamics (4.3) is indeed written in terms of natural sloshing modes (here computed numerically) and the same asymptotic scaling is adopted (see (4.2a–c)).

Nevertheless, the reintroduction of surface tension in the multimodal theory can be challenging. In particular, it is not clear yet how to account for static meniscus and contact line dynamics (Raynovskyy & Timokha 2020). Moreover, in these cases, the element of the Fourier basis (Bessel) functions no more coincide with the actual natural sloshing modes.

The numerical nature of our approach, based on primitive equations, not only allows us to reintroduce straightforwardly surface tension (see Bongarzone et al. 2022a) but also to possibly account (asymptotically) for static meniscus effects and contact angle dynamics while keeping the leading-order dynamics expressed in terms of exact (up to a numerical convergence error) linear natural modes computed numerically. This has been shown possible in a series of works by some of the authors (Viola et al. 2018; Viola & Gallaire 2018; Bongarzone et al. 2021b, 2022b) and makes the present approach in this sense more versatile to study the sloshing problem under non-ideal sidewall conditions.

Most important for the following analysis, is the case of the resonant amplification of higher-order modes. Below a critical liquid depth, typically $H_{cr}\approx 1.05$ for circular cylinders, such amplification (secondary or internal resonances) can happen in the vicinity of the primary resonance. These cases, which require a reordering of the asymptotic scaling, can still be tackled in the framework of the multimodal theory by employing the so-called adaptive modal approach (see chapters 7–9 of Faltinsen & Timokha 2009; Raynovskyy & Timokha 2020). However, secondary resonances may also occur more broadly even far from the primary resonance zone and for $H>H_{cr}$, as in the case of the DC super-harmonic resonance. To the authors’ knowledge, even though its formalization appears possible, no variant of the abovementioned adaptive modal approach capable of dealing with super-harmonic resonances far from the primary one has been reported yet. The WNL analysis of the next section proposes an asymptotic reordering allowing one to deal with the specific case of super-harmonic DC resonances (Bongarzone et al. 2022a).

5. Super-harmonic DC resonance

We now tackle the DC wave response to longitudinal shaking, whose investigation represents the core of the present work. We remind that the DC dynamics occurs at a driving frequency $\varOmega \approx \omega _{2n}/2$ (see figure 4 of Reclari et al. 2014). For the sake of generality, the following analysis is therefore formalized for any mode $(2,n)$, i.e. $\varOmega =\omega _{2n}/2+\lambda$, where $\lambda$ is the small detuning parameter.

By analogy with Bongarzone et al. (2022a), the leading-order solution is here assumed to be given by the sum of a particular solution, given by the linear response to the external forcing, computed by solving (3.5) with $\varOmega =\omega _{2n}/2$ and $m=\pm 1$, and a homogeneous solution, represented by the two natural modes for $(m,n)=(\pm 2,n)$ associated with $\omega _{2n}$, up to their amplitudes to be determined at higher orders (see figure 4). At second order, quadratic terms in $(\varOmega,m)=(\omega _{2n}/2,\pm 1)$ will produce resonant terms in $(\omega _{2n},\pm 2)$. These $\epsilon ^2$-order resonating terms will then require, in the spirit of multiple time-scale analysis, an additional second-order solvability condition, hence suggesting that two slow time scales exist, namely $T_1$ and $T_2$. Thus, the asymptotic scalings of the WNL expansion for DC waves are

(5.1a–d)\begin{equation} f=\epsilon F,\quad \varOmega=\omega_{2n}/2+\epsilon\varLambda, \quad T_1=\epsilon t,\quad T_2=\epsilon^2 t, \end{equation}

with a first-order solution reading

(5.2)\begin{align} \boldsymbol{q}_{1}&=A_2(T_1,T_2)\hat{\boldsymbol{q}}_1^{A_2}\, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}+B_2(T_1,T_2) \hat{\boldsymbol{q}}_1^{B_2}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)} \nonumber\\ &\quad +\frac{1}{2}F\hat{\boldsymbol{q}}_1^{F}\, {\rm e}^{\text{i}((\omega_{2n}/2)t-\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+\frac{1}{2}F\hat{\boldsymbol{q}}_1^{F}\, {\rm e}^{\text{i}((\omega_{2n}/2)t+\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+{\rm c.c.} \end{align}

In (5.2), $\hat {\boldsymbol {q}}_1^{A_2}=\hat {\boldsymbol {q}}_1^{B_2}$, whereas $A_2$ and $B_2$ are the unknown slow time amplitude modulations, here functions of the two time scales $T_1$ and $T_2$. The second-order linearized forced problem reads

(5.3)\begin{align} &(\partial_t\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m}})\boldsymbol{q}_2= \boldsymbol{\mathcal{F}}_2=\boldsymbol{\mathcal{F}}_2^{ij}- \left(\frac{\partial A_2}{\partial T_1}\boldsymbol{\mathsf{B}} \hat{\boldsymbol{q}}_1^{A_2}\, {\rm e}^{\text{i} (\omega_{2n}t-2\theta)} + \frac{\partial B_2}{\partial T_1} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{B_2}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)}+{\rm c.c.}\right) \nonumber\\ &\quad -\text{i}\varLambda F\left( \frac{1}{2} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{F}\, {\rm e}^{\text{i}((\omega_{2n}/2)t-\theta)}{{\rm e}^{\text{i}\varLambda T_1}} +\frac{1}{2} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{F}\, {\rm e}^{\text{i}((\omega_{2n}/2)t+\theta)}{{\rm e}^{\text{i}\varLambda T_1}}+{\rm c.c.}\right). \end{align}

The first-order solution is indeed made of eight different contributions (including the complex conjugates) and it generates, in total, 36 different second-order forcing terms, here implicitly gathered in $\boldsymbol {\mathcal {F}}_2^{ij}$, each characterized by a certain oscillation frequency and azimuthal periodicity. For the sake of brevity, indices $(i,j)$ are used to remind that each forcing is proportional to a quadratic combination of leading-order amplitudes. These indices can assume the following values: $i,j=A_2,B_2,F,\bar {A}_2,\bar {B}_2,\bar {F}$. For instance, the quadratic interaction of $A_2(T_1,T_2)\hat {\boldsymbol {q}}_1^{A_2}\, {\rm e}^{\text {i}(\omega _{2n}t-2\theta )}$ with itself will have indices $(i=A_2,j=A_2)$ and will produce a forcing term proportional to $A_2^2$, i.e. $\boldsymbol {\mathcal {F}}_2^{A_2A_2}$. The additional eight forcing terms, with their complex conjugates, appearing in (5.3), stem from the time derivative of the first-order solution (5.2) with respect to the first-order slow time scale $T_1$. None of the forcing terms in (5.3) are resonant, as their oscillation frequency or azimuthal wavenumber differ from those of the leading-order homogeneous solution, except the two terms produced by the second harmonic of the leading-order particular solution, i.e. $\boldsymbol {\mathcal {F}}_2^{FF}={\frac {1}{4}}F^2 \hat {\boldsymbol {\mathcal {F}}}_2^{FF}\, {\rm e}^{\text {i}(\omega _{2n}t-2\theta )}\, {\rm e}^{\text {i}2\varLambda T_1}+{\frac {1}{4}}F^2\hat {\boldsymbol {\mathcal {F}}}_2^{FF}\, {\rm e}^{\text {i}(\omega _{2n}t+2\theta )}\, {\rm e}^{\text {i}2\varLambda T_1}+{\rm c.c.}$. To avoid secular terms, a second-order compatibility condition is thus imposed, requiring that the following normal form equations are verified:

(5.4a,b)\begin{equation} \frac{\partial A_2}{\partial T_1}=\text{i}\frac{\mu_{{DC}}}{4} F^2\,{\rm e}^{\text{i}2\varLambda T_1},\quad \frac{\partial B_2}{\partial T_1}=\text{i}\frac{\mu_{{DC}}}{4} F^2\,{\rm e}^{\text{i}2\varLambda T_1}. \end{equation}

Taken alone, the dynamics resulting from system (5.4a,b) is still of little relevance, since it can be shown that the wave amplitudes $A_2$ and $B_2$ scale like $\sim {1}/{\varLambda }$, hence diverging symmetrically to infinity for $\varLambda \rightarrow 0$ ($\varOmega \rightarrow \omega _{2n}/2$) in the absence of any restoring term, i.e. the nonlinear mechanism responsible for the finite amplitude saturation, which only comes into play at order $\epsilon ^3$. The expansion must be therefore pursued up to the next order, and thereby one must solve for the second-order solution (Fujimura 1989, 1991).

Figure 4. Spatial structures of the first-order contributions (a) $\boldsymbol {q}_1^F(r,z)\,{\rm e}^{\text {i}}\cos {\theta }$ (SC) and (b) $\hat {\boldsymbol {q}}_1^{A_2}\cos {2\theta }= \hat {\boldsymbol {q}}_1^{B_2}\cos {2\theta }$ (DC) appearing in (5.2) and computed for $t=0$ and $T_1=0$. (c) Superposition of (a) and (b). Here the corresponding amplitudes have been arbitrarily chosen for visualization purposes, but we note that, while amplitudes $A_2$ and $B_2$ still need to be determined, the amplitude of the SC solution (a) is univocally defined once the amplitude, $F$, and the oscillation frequency, $\varOmega$, of the external driving are prescribed.

By substituting (5.2) and (5.4a,b) in the forcing expression, (5.3) can be rewritten as

(5.5)\begin{align} &(\partial_t\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m}})\boldsymbol{q}_2= \boldsymbol{\mathcal{F}}_{2_{NRT}}^{ij}+ \boldsymbol{\mathcal{F}}_{2_{RT}}^{ij}= \boldsymbol{\mathcal{F}}_{2_{NRT}}^{ij}+{\rm c.c.} \nonumber\\ &\quad +\tfrac{1}{4}F^2(\hat{\boldsymbol{\mathcal{F}}}_2^{FF}- \text{i}\mu_{{DC}} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2}) \,{\rm e}^{\text{i}(\omega_{2n}t-2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1}+{\rm c.c.} \nonumber\\ &\quad +\tfrac{1}{4}F^2(\hat{\boldsymbol{\mathcal{F}}}_2^{FF}- \text{i}\mu_{{DC}} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{B_2})\, {\rm e}^{\text{i}(\omega_{2n}t+2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1}+{\rm c.c.}, \end{align}

where the subscripts $_{NRT}$ and $_{RT}$ denote non-resonating and resonating terms, respectively. Note that the term proportional to $\varLambda F$ in (5.3) has been included in the non-resonating forcing terms, while resonant terms are written explicitly. The compatibility condition is now satisfied, meaning that the new resonant forcing term is orthogonal to the adjoint mode, $\hat {\boldsymbol {q}}_1^{A_2{\dagger} }=\bar {\hat {\boldsymbol {q}}}_1^{A_2}$, by construction so that, according to the Fredholm alternative, a non-trivial unique solution can be computed. Hence, we can write the second-order solution as

(5.6)\begin{align} \boldsymbol{q}_2&=(|A_2|^2\hat{\boldsymbol{q}}_2^{A_2\bar{A}_2} +\tfrac{1}{4}|F|^2\hat{\boldsymbol{q}}_2^{F\bar{F}}) \nonumber\\ &\quad +(A_2^2\hat{\boldsymbol{q}}_2^{A_2A_2}\, {\rm e}^{\text{i} (2\omega_{2n}t-4\theta)}+\tfrac{1}{2}\varLambda F\hat{\boldsymbol{q}}_2^{\varLambda F}\, {\rm e}^{\text{i} ((\omega_{2n}/2)t-\theta)}\, {\rm e}^{\text{i}\varLambda T_1} +{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{2}A_2F\hat{\boldsymbol{q}}_2^{A_2F}\, {\rm e}^{\text{i}((3\omega_{2n}/2)t-3\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+\tfrac{1}{2}A_2\bar{F} \hat{\boldsymbol{q}}_2^{A_2\bar{F}}\, {\rm e}^{\text{i} ((\omega_{2n}/2)t-\theta)}\, {\rm e}^{-\text{i}\varLambda T_1}+{\rm c.c.})\nonumber\\ &\quad +(|B_2|^2\hat{\boldsymbol{q}}_2^{B_2\bar{B}_2}+ \tfrac{1}{4} |F|^2\hat{\boldsymbol{q}}_2^{F\bar{F}}) \nonumber\\ &\quad +(B_2^2\hat{\boldsymbol{q}}_2^{B_2B_2}\, {\rm e}^{\text{i} (2\omega_{2n}t+4\theta)}+\tfrac{1}{2}\varLambda F \hat{\boldsymbol{q}}_2^{\varLambda F}\, {\rm e}^{\text{i} ((\omega_{2n}/2)t+\theta)}\, {\rm e}^{\text{i}\varLambda T_1} +{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{2} B_2F\hat{\boldsymbol{q}}_2^{B_2F}\, {\rm e}^{\text{i}((3\omega_{2n}/2)t+3\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+\tfrac{1}{2} B_2\bar{F} \hat{\boldsymbol{q}}_2^{B_2\bar{F}}\, {\rm e}^{\text{i} ((\omega_{2n}/2)t+\theta)}\, {\rm e}^{-\text{i}\varLambda T_1} +{\rm c.c.}) \nonumber\\ &\quad +(A_2B_2\hat{\boldsymbol{q}}_2^{A_2B_2}\, {\rm e}^{\text{i}2\omega_{2n}t}+A_2\bar{B}_2 \hat{\boldsymbol{q}}_2^{A_2\bar{B}_2}\, {\rm e}^{-\text{i} 4\theta}+{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{4}F^2\hat{\boldsymbol{q}}_2^{FF}\, {\rm e}^{\text{i}\omega_{2n}t}\, {\rm e}^{\text{i}2 \varLambda T_1}+\tfrac{1}{4}F\bar{F}\hat{\boldsymbol{q}}_2^{F\bar{F}}\, {\rm e}^{-\text{i}2\theta}+{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{2}A_2F\hat{\boldsymbol{q}}_2^{A_2F}\, {\rm e}^{\text{i}((3\omega_{2n}/2)t-\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+\tfrac{1}{2}A_2\bar{F} \hat{\boldsymbol{q}}_2^{A_2\bar{F}}\, {\rm e}^{\text{i}((\omega_{2n}/2)t-3\theta)}\, {\rm e}^{-\text{i}\varLambda T_1}+{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{2}B_2F\hat{\boldsymbol{q}}_2^{B_2F}\, {\rm e}^{\text{i}((3\omega_{2n}/2)t+\theta)}\, {\rm e}^{\text{i}\varLambda T_1}+\frac{1}{2}B_2\bar{F} \hat{\boldsymbol{q}}_2^{B_2\bar{F}}\, {\rm e}^{\text{i} ((\omega_{2n}/2)t+3\theta)}\, {\rm e}^{-\text{i}\varLambda T_1} +{\rm c.c.}) \nonumber\\ &\quad +(\tfrac{1}{4} F^2\hat{\boldsymbol{q}}_2^{FF}\, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1}+\frac{1}{4} F^2\hat{\boldsymbol{q}}_2^{FF}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1} +{\rm c.c.}). \end{align}

All non-resonant responses in (5.6) are handled similarly, i.e. they are computed in Matlab by performing a simple matrix inversion using standard LU solvers. Although the operator associated with the resonant forcing term, i.e. $(\text {i}\omega _{2n}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2})$, is singular, the value of the normal form coefficient $\mu _{{DC}}$ ensures that a non-trivial solution for $\hat {\boldsymbol {q}}_2^{F^2}$ exists. Diverse approaches can be followed to compute this response, which was here computed by using the pseudo-inverse matrix of the singular operator (Orchini, Rigas & Juniper 2016) obtained via the built-in Matlab function pinv. We also recall that due to the invariant transformation (3.8), only some of the spatial structures appearing in (5.6) need to be computed. Lastly, at third order in $\epsilon$, the problem reads

(5.7)\begin{align} &(\partial_t\boldsymbol{\mathsf{B}}-{\boldsymbol{\mathsf{A}}_{m}})\boldsymbol{q}_3= \boldsymbol{\mathcal{F}}_3 ={-}\frac{\partial A_2}{\partial T_2 }\boldsymbol{\mathsf{B}} \hat{\boldsymbol{q}}_1^{A_2}\, {\rm e}^{\text{i} (\omega_{2n}t-2\theta)}-\frac{\partial B_2}{\partial T_2} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{B_2}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)}\nonumber\\ &\quad - \text{i}\frac{1}{4}2\varLambda F^2\boldsymbol{\mathsf{B}} \hat{\boldsymbol{q}}_2^{F^2}\, {\rm e}^{\text{i} (\omega_{2n}t-2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1} - \text{i}\frac{1}{4}2\varLambda F^2\boldsymbol{\mathsf{B}} \hat{\boldsymbol{q}}_2^{F^2}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)}\, {\rm e}^{\text{i}2\varLambda T_1}\nonumber\\ &\quad + |A_2|^2A_2 \hat{\boldsymbol{\mathcal{F}}}_3^{|A_2|^2A_2}\, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}+ |B_2|^2B_2 \hat{\boldsymbol{\mathcal{F}}}_3^{|B_2|^2B_2}\, {\rm e}^{\text{i}(\omega_{2n}t+2\theta)}\nonumber\\ &\quad + |B_2|^2A_2 \hat{\boldsymbol{\mathcal{F}}}_3^{|B_2|^2A_2}\, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}+ |A_2|^2B_2 \hat{\boldsymbol{\mathcal{F}}}_3^{|A_2|^2B_2}\, {\rm e}^{\text{i}(\omega_{2n}t+2\theta)}\nonumber\\ &\quad +\frac{1}{4}F^2A_2\hat{\boldsymbol{\mathcal{F}}}_3^{|F|^2A_2}\, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}+\frac{1}{4}F^2B_2 \hat{\boldsymbol{\mathcal{F}}}_3^{|F|^2B_2}\, {\rm e}^{\text{i} (\omega_{2n}t+2\theta)}\nonumber\\ &\quad +\frac{1}{4}\varLambda F^2 \hat{\boldsymbol{\mathcal{F}}}_3^{\varLambda F^2} \, {\rm e}^{\text{i}(\omega_{2n}t-2\theta)}\, {\rm e}^{\text{i}2 \varLambda T_1}+\frac{1}{4}\varLambda F^2 \hat{\boldsymbol{\mathcal{F}}}_3^{\varLambda F^2}\, {\rm e}^{\text{i}(\omega_{2n}t+2\theta)}\, {\rm e}^{\text{i}2 \varLambda T_1}+\text{N.R.T.}+{\rm c.c.}, \end{align}

where the first two forcing terms arise from the time derivative of the first-order solution with respect to the second-order slow time scale $T_2$ and from that of the second-order solution with respect to the first-order slow time scale $T_1$, respectively (see Appendix D of Bongarzone et al. (2022a) for the full expression of $\boldsymbol {\mathcal {F}}_2$ and $\boldsymbol {\mathcal {F}}_3$). Once again, all terms explicitly written in (5.7) are resonant, as they share the same pair $(\omega _{2n},\pm 2)$ as the first-order homogeneous solutions, hence a third-order compatibility condition, leading to the following normal form, must be enforced,

(5.8a)\begin{gather} \frac{\partial A_2}{\partial T_2}=\text{i}\frac{\zeta_{{DC}}}{4} \varLambda F^2\,{\rm e}^{\text{i}2\varLambda T_1}+\text{i} \frac{\chi_{{DC}}}{4}A_2F^2+\text{i}\nu_{{DC}} |A_2|^2A_2+ \text{i}\xi_{{DC}}|B_2|^2A_2, \end{gather}

(5.8b)\begin{gather}\frac{\partial B_2}{\partial T_2}=\text{i}\frac{\zeta_{{DC}}}{4} \varLambda F^2\,{\rm e}^{\text{i}2\varLambda T_1}+\text{i}\frac{\chi_{{DC}}}{4}B_2F^2+\text{i}\nu_{{DC}} |B_2|^2B_2+\text{i}\xi_{{DC}}|A_2|^2B_2. \end{gather}

where the coefficients are defined in Appendix A.

As a last step in the derivation of the final amplitude equation for the DC waves and in order to eliminate the implicit small parameter $\epsilon$, we unify systems (5.4a,b), (5.8a) and (5.8b) into a single system of equations recast in terms of the physical time $t=T_1/\epsilon =T_2/\epsilon ^2$, physical forcing control parameters, $f=\epsilon F$, $\lambda =\epsilon \varLambda$ and total amplitudes, $A=\epsilon A_2\, {\rm e}^{-\text {i}2\lambda t}$ and $B=\epsilon B_2\, {\rm e}^{-\text {i}2\lambda t}$. This is achieved by summing (5.4a,b)–(5.8a) and (5.8b) along with their respective weights $\epsilon ^2$ and $\epsilon ^3$, thus obtaining

(5.9a)\begin{gather} \frac{{\rm d} A}{{\rm d} t}={-}\text{i}\left(2\lambda-\frac{\chi_{{DC}}}{4} f^2\right) A + \text{i}\frac{\left(\zeta_{{DC}}\lambda + \mu_{{DC}}\right)}{4} f^2 + \text{i}\nu_{{DC}} |A|^2A+\text{i}\xi_{{DC}}|B|^2A, \end{gather}

(5.9b)\begin{gather}\frac{{\rm d} B}{{\rm d} t}={-}\text{i}\left(2\lambda-\frac{\chi_{{DC}}}{4} f^2\right) B + \text{i}\frac{\left(\zeta_{{DC}}\lambda + \mu_{{DC}}\right)}{4} f^2 + \text{i}\nu_{{DC}} |B|^2B+\text{i}\xi_{{DC}}|A|^2B. \end{gather}

We note that no second-order homogeneous solutions, e.g. proportional to amplitudes $C_2(T_1,T_2)$ and $D_2(T_1,T_2)$, have been accounted for in (5.6), as their presence will produce two resonant third-order terms, $({\partial C_2}/{\partial T_1})\boldsymbol{\mathsf{B}}\hat {\boldsymbol {q}}_2^{C_2}\, {\rm e}^{\text {i}(\omega _{2n}t-2\theta )}$ ($\hat {\boldsymbol {q}}_2^{C_2}=\hat {\boldsymbol {q}}_2^{A_2}$) and $({\partial D_2}/{\partial T_1})\boldsymbol{\mathsf{B}}\hat {\boldsymbol {q}}_2^{D_2}\, {\rm e}^{\text {i}(\omega _{2n}t+2\theta )}$ ($\hat {\boldsymbol {q}}_2^{C_2}=\hat {\boldsymbol {q}}_2^{A_2}$), that can be incorporated in the final amplitude equations (5.9a) and (5.9b) by simply defining $A=\epsilon ( A_2+\epsilon C_2)\, {\rm e}^{-\text {i}2\lambda t}$ and $B=\epsilon ( B_2+\epsilon D_2)\, {\rm e}^{-\text {i}2\lambda t}$.

As in § 4, we first turn to polar coordinates, $A=|A|{\rm e}^{\text {i}\varPhi _A}$ and $B=|B|{\rm e}^{\text {i}\varPhi _A}$, and we split the modulus and phase parts of (5.9a) and (5.9b). We then look for stationary solutions, $d/dt=0$ with $|A|,\,|B|\ne 0$ ($\varPhi _A=\varPhi _B=\varPhi =0,{\rm \pi}$; see § 4). By summing and subtracting (5.9a) and (5.9b), after introducing the auxiliary amplitudes $|a|=|A|+|B|$ and $|b|=|A|-|B|$, the following implicit relations are obtained:

(5.10a)\begin{gather} f^2=|a|\left(2\lambda-\frac{\nu_{{DC}}+\xi_{{DC}}}{4}|a|^2- \frac{3\nu_{{DC}}-\xi_{{DC}}}{4}|b|^2\right) \frac{4}{(|a|\chi_{{DC}}\pm2(\zeta_{{DC}}\lambda+ \mu_{{DC}}))}, \end{gather}

(5.10b)\begin{gather}0=|b|\left(\frac{\chi_{{DC}}}{4}f^2-\left(2\lambda-\frac{\nu_{{DC}}+ \xi_{{DC}}}{4}|b|^2-\frac{3\nu_{{DC}}-\xi_{{DC}}}{4}|a|^2 \right)\right). \end{gather}

Here $f=a_x\varOmega ^2$ and $\lambda =\varOmega -\omega _{2n}/2$. By analogy with harmonic forcing conditions, two possible (super-harmonic) solutions exist, i.e. a planar wave solution for $|b|=0$,

(5.11)\begin{equation} f=\sqrt{|a|\left(2\lambda-\frac{\nu_{{DC}}+\xi_{{DC}}}{4}|a|^2\right) \frac{4}{(|a|\chi_{{DC}}\pm2(\zeta_{{DC}}\lambda+\mu_{{DC}}))}}, \end{equation}

and a swirling solution for $|b|\ne 0$, defined by

(5.12a)\begin{gather} |b|^2=\left(2\lambda-\frac{\chi_{{DC}}}{4}f^2- \frac{3\nu_{{DC}}-\xi_{{DC}}}{4}|a|^2\right) \left(\frac{4}{\nu_{{DC}}+\xi_{{DC}}}\right), \end{gather}

(5.12b)\begin{gather}f=\sqrt{2|a|\left(\frac{\xi_{{DC}}-\nu_{{DC}}}{\nu_{{DC}}+ \xi_{{DC}}}\right)(2\lambda-\nu_{{DC}}|a|^2) \frac{4}{\left(2|a|\dfrac{(\xi_{{DC}}-\nu_{{DC}} )}{(\nu_{{DC}}+\xi_{{DC}})}\chi_{{DC}} \pm2(\zeta_{{DC}}\lambda+\mu_{{DC}})\right)}}, \end{gather}

where only real solutions corresponding to $f=a_x\varOmega ^2>0$ are retained, as the combinations $a_x\varOmega ^2<0$ are not physically meaningful.

The stability of such stationary solutions $\boldsymbol {y}_{s}=(|A|,\varPhi _A,|B|,\varPhi _B)$ is computed by introducing small amplitude and phase perturbations ($\ll 1$) with the ansatz $\boldsymbol {y}_{p}(t)=(|A_p|,\varPhi _{A,p},|B_p|,\varPhi _{B,p})\, {\rm e}^{st}$ in (5.9a) and (5.9b), which are then linearized around $\boldsymbol {y}_0$, hence obtaining at first order an eigenvalue problem in the complex eigenvalue $s=s_R+\text {i}s_I$. For each $(|A|,\varPhi _A,|B|,\varPhi _B)$, one obtains four eigenvalues $s$ and if the real part $s_R$ of at least one of these eigenvalues is positive, then that configuration is deemed as unstable. An analogous procedure has been followed for the case of harmonic resonances discussed in § 4.

Once the various branches for $|a|$ and $|b|$ as a function of $\tau =\varOmega /\omega _{2n}$ and at a fixed non-dimensional shaking amplitude $a_x$ are computed and their stability is determined, amplitudes $A$ and $B$ are substituted in (5.2) and (5.6), so that the total flow solution predicted by the WNL for DC waves is reconstructed as

(5.13)\begin{equation} \boldsymbol{q}_{DC}=\{\varPhi,\eta\}^T=\epsilon\boldsymbol{q}_1+ \epsilon^2\boldsymbol{q}_2. \end{equation}

As discussed in Bongarzone et al. (2022a) for circular sloshing, although the quantitative dependence on the external control parameters, i.e. driving amplitude and frequency, is different with respect to the SC case, e.g. $f^2$ instead of $f$, system (5.9a) and (5.9b) is essentially analogous to that given in (4.7a) and (4.7b). Indeed, (5.9a) and (5.9b) contain four main contributions,

(5.14a–d)\begin{equation} \lambda \leftrightarrow \left(2{\lambda}-\frac{\chi_{{DC}}}{4}f^2\right),\quad \mu_{{SC}}f\leftrightarrow \frac{\zeta_{{DC}}\lambda+\mu_{{DC}}}{4}f^2,\quad \nu_{{SC}}\leftrightarrow\nu_{{DC}},\quad \xi_{{SC}}\leftrightarrow \xi_{{DC}}, \end{equation}

corresponding respectively to a detuning term (forcing amplitude dependent), an additive (quadratic) forcing term (driving frequency dependent), the classic cubic restoring term and, lastly, the cubic term dictating the nonlinear interaction between the two counter-propagating travelling waves. For these reasons, figure 5 shows the nonlinear amplitude saturation for $|a|=|A|+|B|$ and $|b|=|A|-|B|$ that are reminiscent of those commented and displayed by Faltinsen et al. (2016) in their figure 7 with regard to harmonic system responses, although the phases associated to each super-harmonic branch are ${\rm \pi}$ shifted with respect to the their harmonic analogous.

Figure 5. Typical response curve for $a$ and $b$ for a fluid depth $H=1.5$ with longitudinal super-harmonic forcing of amplitude $a_x=0.2$. Panel (a) shows a projection of the three-dimensional branch structure $(\varOmega /\omega _{21},|a|,|b|)$ in the $(\varOmega /\omega _{21},|a|)$ plane, whereas panel (b) shows the same projection, but on the $(\varOmega /\omega _{21},|b|)$ plane. Black solid lines mark stable steady-state planar waves, whereas light blue solid lines indicate stable steady-state swirling waves. Dashed lines denote the corresponding unstable branches. Here U denotes the turning point, H denotes the Hopf bifurcation and P denotes the Poincaré bifurcation. For completeness, the phase values $\varPhi _A=\varPhi _B=\varPhi =0$ or ${\rm \pi}$ associated to each branch are reported in panel (a).

A more detailed description of the bifurcation diagrams shown in figures 5(a) and 5(b) is given in Faltinsen et al. (2016). Here we limit to note that the branching diagrams contain three bifurcation points, namely U (turning point), H (Hopf bifurcation) and P (Poincaré bifurcation, Miles 1984b), whose positions determine the frequency ranges where stable planar (standing), swirling or irregular waves are theoretically expected. By keeping track of the position of these three bifurcation points in the $(\varOmega /\omega _{21},|a|)$ plane as the forcing amplitude, $a_x$, is varied, one can draw a super-harmonic stability chart in the $(\varOmega /\omega _{21},a_x)$ plane similar to that of figure 3 for harmonic resonances and which is shown in figure 6.

Figure 6. Estimates of bounds, in the $(\varOmega /\omega _{21},a_x)$ plane, between the frequency ranges where planar, irregular and swirling waves occur when the container undergoes a longitudinal and super-harmonic motion at a forcing frequency $\varOmega \approx \omega _{21}/2$. In this range of frequency, the theory predicts the superposition of an unconditionally stable planar SC wave ($m=\pm 1$) oscillating harmonically with the driving frequency and a super-harmonic DC dynamics ($m=\pm 2$), which can manifest itself via planar, swirling or irregular wave motions. The stability boundaries (black solid lines) were computed for a fluid depth $H=1.5$, as in Royon-Lebeaud et al. (2007). The corresponding values of the normal form coefficients appearing in (5.9a) and (5.9b) are given in table 1.

The first striking difference with respect to the harmonic stability chart of figure 2 is the opposite curvature of the stability boundaries between the various super-harmonic regimes. As mentioned above, this is due to the quantitative dependence of the additive forcing term in system (5.9a) and (5.9b) on the driving amplitude, which is here quadratic in $f$, thus leading to the square root in (5.11) (planar DC) and (5.12b) (swirling DC).

Furthermore, there is a substantial difference in terms of free-surface patterns. As suggested by the form of the first-order solution (5.2), the leading-order dynamics, governing the super-harmonic system response to longitudinal forcing, results from a superposition of a stable planar (or standing) SC wave, oscillating harmonically at a frequency $\omega _{{SC}}=\varOmega \approx \omega _{2n}/2$ and generated by the two $m=\pm 1$ counter-rotating travelling waves of equal amplitudes, and a super-harmonic DC wave dynamics oscillating at a frequency of approximately $\omega _{{DC}}=2\varOmega \approx \omega _{2n}$ (period halving). We can also relate the SC wave frequency to its own natural frequency by writing $\omega _{{SC}}=\varOmega \approx \frac {1}{2}\sqrt {{k_{2n}\tanh {k_{2n}H}}/{k_{1n}\tanh {k_{1n}H}}}\omega _{1n}$, which is $\approx 0.657\omega _{1n}$ (in deep water) for $n=1$ and approaches $0.5\omega _{1n}$ for large $n$, hence showing that the DC resonance always occurs far from the primary harmonic resonance.

When the amplitudes of the two travelling waves with $m=\pm 2$ are equal, i.e. $|A|=|B|$ (or $|b|=0$), the DC dynamics manifests itself via planar motion and the global solution takes the form of a planar wave (planar ${\rm SC}+{\rm DC}$, light blue shaded region in figure 6). On the contrary, when $|A|\ne |B|\ne 0$, one of the two $m=\pm 2$ waves dominates over the other and a stable swirling motion, responsible for the system symmetry breaking, is established. In this case, the total solution is given by the sum of a harmonic planar SC wave and a super-harmonic swirling DC wave (swirling DC + planar SC, green shaded region in figure 6). The white dotted region and the light red shaded regions in figure 6 correspond, respectively, to the super-harmonic irregular motion regime (see § 6 for further details) and to the multi-solution range where both types of motion are possible depending on the initial conditions, i.e. to the region of hysteresis.

6. Experiments

In this section we present our experimental set-up dedicated to the generation and characterization of sloshing waves under longitudinal super-harmonic forcing with driving (dimensionless) frequency $\varOmega \approx \omega _{21}/2$. The bounds between the different regimes for the resulting super-harmonic wave are experimentally retrieved as a function of the driving amplitude and frequency, and compared with the theoretical estimates. Finally, we measure the wave amplitude saturation in the vicinity of the super-harmonic resonance, and compare it with the theoretical WNL prediction (5.13).

6.1. Experimental set-up

The experimental set-up used to generate the sloshing waves in the cylindrical container and to observe the resulting free-surface motion is shown in figure 7. A Plexiglas cylindrical container of height 50 cm and inner diameter $D=2R= 17.2$ cm, partially filled with a column of distilled water of height $h= 11$ cm, is fixed on a single-axis linear motion actuator (AEROTECH PRO165LM). Sloshing waves are generated by imposing to the container a longitudinal sinusoidal forcing of angular frequency $\bar {\varOmega }$ and amplitude $\bar {a}_x$. The motion of the fluid free surface is recorded with a digital camera (NIKON D850) coupled with a Nikon 60 mm f/2.8D lens and operated in slow motion mode, allowing for an acquisition frequency of 120 frames per second. The optical axis of the camera is aligned with the container motion axis. A LED panel (not depicted in figure 7) placed behind the tank provides back illumination of the fluid free surface for a better optical contrast.

Figure 7. Experimental apparatus.

The actuation of the moving stage as well as the camera triggering for movie recording are set and controlled via a home-made Labview program. In a typical experiment, the container undergoes a harmonic motion of fixed amplitude in the range $4\ {\rm mm} \leq \bar {a}_x \leq {\rm 34}\ {\rm mm}$ (i.e. $a_x = \bar {a}_x /R \in [0.05, 0.40]$), while a sweep in forcing frequency is implemented within the interval $\bar {\varOmega }/ 2{\rm \pi} \in [1.35\ {\rm Hz}, 1.58\ {\rm Hz}]$ corresponding to the dimensionless range $\varOmega / \omega _{21} \in [0.45, 0.53]$. Each frequency step lasts 100 oscillation periods while the frequency increment between two consecutive steps is typically of 10 mHz. Along the sweeping, a movie is recorded for each $(\bar {a}_x, \bar {\varOmega })$ set of parameters. To ensure that the steady-state amplitude regime is established at each step in the recorded free-surface dynamics, the camera is triggered only after a certain number of cycles, typically 50, see Appendix C.

6.2. Analysis of the free-surface dynamics

6.2.1. Qualitative observations

While operating a sweep in forcing frequency at a fixed forcing amplitude, we observe in the vicinity of the super-harmonic resonance three different kinds of motion, namely planar, irregular and swirling ones, whose occurrence depends on the forcing amplitude and frequency; see, for instance, the snapshots displayed on figure 8 or the videos provided among the supplementary movies available at https://doi.org/10.1017/jfm.2023.438.

Figure 8. Images of the fluid free surface while the container is subjected to a longitudinal harmonic forcing of amplitude $a_x=\bar {a}_x/R \approx 0.23$ at various driving angular frequencies $\varOmega$ close to $\omega _{21}/2$. The fluid free surface is observed in the direction aligned with the container motion. For each driving frequency (a), (b) and (c), the time interval between two snapshots is about $T/4$, with $T=2 {\rm \pi}/ \varOmega$ the corresponding oscillation period. On each snapshot, the vertical middle axis is represented by a red dotted line. For a forcing frequency $\varOmega \approx 0.48 \omega _{21}$ (a) and $\varOmega \approx 0.52 \omega _{21}$ (c) the free-surface image at each time $t$ is mirror symmetric with respect to the middle vertical axis, the signature of a planar wave regime, while the symmetry is broken for $\varOmega \approx 0.50 \omega _{21}$ (c) revealing a swirling wave. Results are shown for (a) $\varOmega / \omega _{21}\approx 0.48$, (b) $\varOmega / \omega _{21}\approx 0.50$ and (c) $\varOmega / \omega _{21}\approx 0.52$.

For a given (and large enough) amplitude and starting from a frequency higher than a certain amplitude-dependent threshold $\varOmega _P(a_x)$, the free surface responds to the longitudinal harmonic forcing by displaying a planar dynamics such as shown in figure 8(c). When the critical frequency $\varOmega =\varOmega _P(a_x)$ is reached, the motion bifurcates to a swirling wave, which propagates along the container wall with a stationary amplitude; see figure 8(b). The wave can rotate either clockwise or anti-clockwise (both rotation directions were observed along the experiments). When the forcing frequency is further decreased below a critical frequency $\varOmega = \varOmega _H(a_x) \approx \omega _{21}/2 < \varOmega _P(a_x)$, the free surface exhibits an irregular dynamics, characterized by a switching between planar and swirling motion (not shown in figure 8). For forcing frequencies lower than a certain threshold $\varOmega < \varOmega _U(a_x)$, the free-surface motion stabilizes into a steady planar wave such as shown in figure 8(a). All together, these observations are qualitatively consistent with the outcomes of the WNL analysis in § 5, that predicts the existence of three different dynamical regimes – namely planar, irregular and swirling motions – for a longitudinal forcing frequency in the vicinity of $\omega _{21}/2$. One of the main purposes of the present experimental investigation is to determine the amplitude-dependent frequency bounds of these different regimes and to compare them to our theoretical prediction of the positions of the bifurcation points U (turning point), H (Hopf bifurcation) and P (Poincaré bifurcation) (see figures 5 and 6).

6.2.2. Procedure

Since the camera optical axis is aligned with the direction of the container motion, we note that a planar wave is characterized by its symmetry with respect to the vertical middle axis of the container image, whereas a swirling wave breaks this symmetry while travelling clockwise or anti-clockwise along the container walls; see figure 8.

We take benefit of these observations to build a more quantitative description of the free-surface dynamics, with the aim of identifying the various types of sloshing waves in the vicinity of the super-harmonic resonance. This can be done by exploiting the symmetry properties of the image of the free-surface response with respect to the vertical middle axis of the container image, and by characterizing the regularity of these waves as a function of the forcing parameters, so as to identify the irregular regime.

To do so, the time evolution of the free-surface dynamics is extracted from the movies along vertical directions that are mirror symmetric with respect to the vertical middle axis of the container image. Comparing the resulting temporal signals with each other allows one to discriminate between planar and swirling motions and to study the wave regularity.

The first step is to attach to each frame $t$ of a given movie, a Cartesian reference frame $(Y(t), Z(t))$, such that $Y(t)=0$ corresponds to the vertical middle axis of the container image, and that $Y(t)=R$ represents the right-hand side edge of the container image. To this end, the edges of the container are automatically detected in a dedicated Matlab program. The vertical $Z(t)$ axis ($Y(t)=0$) on the frame corresponding to time $t$ is then set as the middle line between these two edges, while the distance between both edges sets the scale of the horizontal direction $Y$. Note that we neglect the 4 mm thickness of the container wall.

A direction $y \in [-R, R]$ is then chosen to extract from each frame corresponding to time $t_i$, the intensity profile $I_{t_i}(y)$ along the vertical line $Y(t_i)=y$. The resulting intensity profiles are then plotted as a function of time to build an image $I(y)$ composed as $I(y)=[I_{t_1}(y), I_{t_2}(y),\ldots ]$, such as displayed in figure 9(c).

Figure 9. General procedure for the analysis of the free-surface dynamics. (a) On each frame, the edges of the container are detected (black dotted lines) and the vertical $Z(t)$ axis is set as the middle line between these two edges, while the scale of the horizontal direction $Y(t)$ is fixed by the distance between both edges. (b) Schematic of the container illustrating the link between the Cartesian coordinate system $(Y(t), Z(t))$ attached to each frame and the cylindrical coordinates in the referential frame of the container. (c) Left: the intensity profile along a vertical line of coordinate $(Y(t)=y)$ with $y \in [-R, R]$ is then measured on each frame $t$ and plot as a function of time (here for $y=0$). The position of the front contact line at the azimuthal coordinate $\theta =\arcsin (y/R)=0$ as a function of time is highlighted in red. Right: frames from which the intensity profiles at times $t_i$ and $t_j$ on the left-hand side image, along the line $(Y(t)=0)$ (represented by a red dotted line), are extracted. At $t_i$, the wave is climbing the front wall of the container (with respect to the camera position) whereas at $t_j$ it reaches the back of the tank.

We note that at each time $t$, the intensity profile $I_{t}(y)$ contains the intersection of the front contact line image with the vertical axis $(Y(t)=y)$, that corresponds to the point of coordinates $(R, \theta, \eta (R, \theta, t))$ in the moving cylindrical frame of reference of the container, where $\theta =\arcsin (y/R)$ (see figure 9b). As a consequence, the final image $I(y)$ also contains the dynamics of the front contact line in the azimuthal direction $\theta$.

The resulting image $I(y)$ exhibits a periodic dark pattern that represents the free-surface response to the harmonic forcing, see an example in figure 9(c) in which $y=0$. Indeed, on each frame of the movie, the free surface appears as the darkest feature, so that the intensity profile along a given line $(Y(t)=y)$ actually represents the vertical extension of the free surface at time $t$ along this direction, which is maximal whenever the sloshing wave reaches its maximal elevation $\max _t \eta ( R, \theta, t)$ along the azimuthal direction $\theta = \arcsin (y/R)$ (in the front of the container with respect to the camera position, corresponding to $\theta \in ]-{\rm \pi} /2, {\rm \pi}/2[$) or along $\theta = {\rm \pi}- \arcsin (y/R)$ (in the back of the container). Furthermore, when the contact line reaches its maximal elevation in the front of the container, the free surface is imaged from below, so that it appears darker than when the maximal elevation is reached in the back, where the free surface is imaged from above; see the snapshots in figure 9(c). These observations allow us to identify in the image $I(y)$ the position as a function of time of the front contact line $\eta ( R, \theta, t)$, with $\theta =\arcsin (y/R)$, as highlighted in red in figure 9(c), and following the method detailed in Appendix D.

Note that this procedure does not give a quantitative access to the actual amplitude of the front contact line oscillations, since the intensity profiles $I_{t_i}(y)$ constituting the image $I(y)$ are simply juxtaposed with each other without rescaling the pixel width along the vertical direction. However, the position extracted from $I(y)$ of the image of the points of coordinates $\eta (R, \pm \theta, t)$ as a function of time still encloses the symmetry properties of the free-surface response, its regularity as well as its frequency content, which are the only quantities needed in order to identify the wave regimes.

6.3. Regularity and frequency content of the free-surface response

The resulting image $I(y)$ is then revealing of the free-surface dynamics $\eta (r, \theta, t)$ and, in particular, of its dynamics at the front wall $\eta (r=R, \theta =\arcsin (y/R), t)$. Figure 10(a–d) displays $I(0)$ for various forcing frequencies close to the super-harmonic resonance, at the same forcing amplitude. These images reveal that depending on the forcing frequency, the free-surface oscillations (dark periodic pattern) can be either regular (a), (c) and (d) – i.e. the oscillations are enclosed into an envelope of constant amplitude – or irregular (b) with a temporal modulation of the amplitude envelope. Therefore, the profiles $I(0)$ allow us to characterize the regularity of the sloshing wave and, in particular, to identify the irregular regime. The latter will be described in more detail in § 6.6, but such details are not needed for the identification of the irregular regime bounds, for which the analysis of the regularity property of the $I(0)$ pattern is sufficient. Therefore, in the following we will focus on the regular planar and swirling motions, which cannot be unambiguously distinguished from each other on the basis of the profiles $I(0)$.

Figure 10. (a–d) Intensity profiles as a function of time along the vertical middle axis $(Y=0)$ – denoted $I(0)$ – for various forcing frequencies $\varOmega$ in the vicinity of the super-harmonic resonance $\varOmega \approx \omega _{21}/2$ at the same forcing amplitude $\bar {a}_x/R \approx 0.23$. In each case, the profile $I(0)$ is extracted from a movie whose recording has been started after about 50 oscillation cycles following each change in forcing frequency, thus ensuring that initial transients are filtered out (see also Appendix C). (e–h) Power spectral densities (PSD) – normalized by the maximal peak amplitude – corresponding to the front contact line dynamics as extracted from the profiles $I(0)$ displayed in panels (a)–(d).

Figure 10(e–h) displays the (normalized) power spectral densities of the front contact line dynamics $\eta (R, \theta =0, t)$ extracted from the profiles $I(0)$ (a–d). It appears that in all cases, the energy of the sloshing wave is massively distributed to its first (harmonic) and second (super-harmonic) component, while the contribution of higher modes is fairly negligible. This incidentally implies that the symmetry properties of a regular wave are directly linked to the symmetry properties of these two first oscillation modes.

In other words, a planar dynamics should necessarily consist in the superposition of two planar waves: a planar SC wave harmonically oscillating with the driving frequency $\varOmega$ and one super-harmonic planar DC wave oscillating at $\omega _{DC}= 2 \varOmega \approx \omega _{21}$. On the other hand, a swirling dynamics must contain at least one symmetry-breaking (swirling) component that, as predicted by the present WNL analysis, should correspond to the super-harmonic $\omega _{21}$ component.

6.4. Symmetry properties of the regular regimes: planar versus swirling waves

We now focus on the regular regimes, namely the steady planar and swirling motions. As stated before, the profiles $I(0)$ cannot discriminate between a planar and a swirling dynamics and instead only contain information on their regularity and their frequency content. To distinguish a planar from a swirling motion, we then compare the profiles along two ($Y \neq 0$) directions that are symmetric with respect to the vertical middle axis of the container image.

Figures 11(b), 11(d) and 11(f) show composite images, each produced using the Matlab function imshowpair applied to the pair $I(R/2)$ and $I(-R/2)$ for three different forcing frequencies that both result in a regular motion (the same forcing parameters as in figure 10a,c,d). Briefly, imshowpair $(I_\alpha,I_\beta )$ creates from a pair of greyscale images $I_\alpha$ and $I_\beta$, a red, green and blue (RGB) image where each pixel is represented by a RGB triplet, the R intensity being the intensity of the corresponding pixel in $I_\alpha$, and the G and B intensities being equal to the intensity of the corresponding pixel in $I_\beta$. A pixel where $I_\alpha$ and $I_\beta$ have the same intensity will be represented by a RGB triplet of the form $[a, a, a]$, where $a \in [0, 255]$, i.e. will appear as grey. On the contrary, if this pixel has a much larger intensity on $I_\alpha$ (respectively on $I_\beta$) than it has on $I_\beta$ (respectively on $I_\alpha$), it will appear in red (respectively in cyan) on the resulting composite image. The composite images displayed in figure 11(b,d,f) thus highlight in each case the differences between $I(R/2)$ and $I(-R/2)$. They are then a direct signature of the symmetry of the free-surface dynamics with respect to the vertical middle axis $(Y=0)$, and reveal two different kinds of motion: (i) a planar motion, for which $I(R/2)$ and $I(-R/2)$ perfectly overlap with each other due to the mirror symmetry of the wave; and (ii) a circular motion, characterized by a symmetry breaking between the right- and left-hand side free-surface dynamics: the maximum of the wave along $\theta =\arcsin (1/2)={\rm \pi} /6$ is indeed phase shifted with respect to the maximum of the wave along $\theta =-{\rm \pi} /6$, thus revealing a travelling wave propagating along the wall of the container.

Figure 11. Symmetry properties of the stationary waves. (a,c,e) Images of the fluid free surface while the container is subjected to a longitudinal harmonic forcing of amplitude $\bar {a}_x/R \approx 0.23$ at various driving angular frequencies $\varOmega$ close to $\omega _{21}/2$: (a) $\varOmega \approx 0.48 \omega _{21}$, (c) $\varOmega \approx 0.50 \omega _{21}$ and (e) $\varOmega \approx 0.52 \omega _{21}$ (same forcing parameters as in figure 10a,c,d). For each driving frequency (a,c,e), the time interval between two snapshots is about $T/4$, with $T=2 {\rm \pi}/ \varOmega$ the corresponding oscillation period. On each snapshot, the vertical axes $(Y=R/2)$ and $(Y=-R/2)$ are represented by a blue and red dotted line, respectively. For a forcing frequency $\varOmega =0.48 \omega _{21}$ (a) and $\varOmega =0.52 \omega _{21}$ (e) the free-surface image at each time $t$ is mirror symmetric with respect to the middle vertical axis, while the symmetry is broken for $\varOmega = 0.50 \omega _{21}$. (b,d,f) Superposition of the intensity profiles as a function of time along the vertical axis $(Y=R/2)$ and $(Y=-R/2)$ – denoted $I(R/2)$ (in blue) and $I(-R/2)$ (in red), respectively – for the same forcing parameters as in figure 10(a, c, e). The grey regions show where $I(R/2)$ and $I(-R/2)$ have the same intensities.

To determine which $\omega$ component is responsible for the symmetry breaking induced by the swirling motion, we extract from $I(R/2)$ and $I(-R/2)$ the position of the front contact line as a function of time $\eta (R, \theta, t)$, where $\theta = \pm {\rm \pi}/6$; see figure 12(b,e). This makes it possible to compute the power spectrum of both signals, as well as the phase difference between the phase angle of their components that oscillate at the frequencies corresponding to their spectrum's first and second peaks (see figure 12c–f). A planar wave oscillating at a frequency $\omega$ is then characterized by the $\omega$ components of $\eta (R, {\rm \pi}/6, t)$ and of $\eta (R, -{\rm \pi} /6, t)$ being in phase with each other, while a swirling wave is characterized by a $m {\rm \pi}/3$-phase shift between the $\omega$ components of these signals, where $m$ denotes the azimuthal wavenumber of the swirling wave ($m=1$ for an harmonically oscillating SC wave, $m=2$ for a super-harmonic DC wave).

Figure 12. Analysis of the steady-state free-surface dynamics under an harmonic forcing of amplitude $a_x=0.23$ and frequency (a–c) $\varOmega /\omega _{21} \approx 0.48$ and (d–f) $\varOmega / \omega _{21} \approx 0.50$. (a,d) Image of the free surface with vertical lines intersecting the image of the front contact line at the points of coordinates $(R, {\rm \pi}/6, \eta (R, {\rm \pi}/6, t))$ (blue arrows) and $(R, -{\rm \pi} /6, \eta (R, -{\rm \pi} /6, t))$ (red arrows) in the moving reference frame of the container. (b,e) Normalized elevation of the front contact line $\tilde {\eta }(R, {\rm \pi}/6, t)$ (blue dots) and $\tilde {\eta }(R, -{\rm \pi} /6, t)$ (red dots) extracted from the corresponding profiles $I(R/2)$ and $I(-R/2)$ (not shown here). The $\tilde {\eta }$ functions are defined according to $\tilde {\eta }(R, \theta, t)= (\eta (R, \theta, t)-\sigma ) / \delta$, where $\sigma =\min _t(\eta (R, \theta, t))+\max _t(\eta (R, \theta, t) )/2$ and $\delta =\max _t(\eta (R, \theta, t))-\min _t(\eta (R, \theta, t))$. (c,f) Left: power spectral densities of $\tilde {\eta }(R, {\rm \pi}/6, t)$ (blue dots) and of $\tilde {\eta }(R, -{\rm \pi} /6, t)$ (red dots). Right: absolute value of the phase shift between the components of $\tilde {\eta }(R, {\rm \pi}/6, t)$ and of $\tilde {\eta }(R, -{\rm \pi} /6, t)$ oscillating at the frequencies corresponding to the first peak ($\omega = \varOmega$) and to the second peak ($\omega = 2 \varOmega$) of the power spectra.

The Fourier analysis of the signals $\eta (R, {\rm \pi}/6, t)$ and $\eta (R, -{\rm \pi} /6, t)$ reveals that for forcing frequencies $\varOmega$ close to $\omega _{21}/2$, the free-surface motion mostly results from the combination of a SC wave harmonically oscillating at the forcing frequency $\omega _{SC}=\varOmega \approx \omega _{21}/2$, and of a super-harmonic DC wave oscillating at a frequency $\omega _{DC}=2 \varOmega \approx \omega _{21}$. From figures 12(c) and 12(f), it is clear that the SC wave is a planar wave for both planar (figure 12c) and swirling (figure 12f) dynamics, as revealed by the vanishing phase shift between the harmonic components of $\eta (R, {\rm \pi}/6, t)$ and $\eta (R, -{\rm \pi} /6, t)$ in both cases. On the other hand, the phase shift between the super-harmonic components is zero in the case of the planar dynamics, and close to $2 {\rm \pi}/3$ in the case of the swirling dynamics, the signature of a DC swirling wave. These observations are general to the whole range of forcing frequencies and amplitudes investigated along this study: in the vicinity of the super-harmonic resonance, the SC wave is always a planar wave, as revealed by the vanishing phase shift between the harmonic components of $\eta (R, {\rm \pi}/6, t)$ and $\eta (R, -{\rm \pi} /6, t)$ for both planar and swirling dynamics (this is also true in the case of the irregular regime, see later § 6.6). In the case of a regular dynamics, the DC wave is either a planar (for vanishing phase shift between the corresponding components) or a swirling wave (characterized by a $2 {\rm \pi}/3$ phase shift between the $\omega _{21}$ component of the right- and left-hand side signals), depending on the exact ratio between $\varOmega$ and $\omega _{21}$, as well as on the forcing amplitude $\bar {a}_x/R$.

6.5. Experimental estimate of regime bounds

From the above analysis, it appears that consistently with the predictions provided by our theoretical WNL analysis, the sloshing waves resulting from the longitudinal super-harmonic forcing of the container at a frequency $\varOmega \approx \omega _{21}$, consist in the superposition of a planar SC wave, harmonically oscillating with the forcing at $\omega =\varOmega$, and of a DC wave, which can exhibit either a planar, irregular or swirling dynamics.

Having identified the three different regimes for the free-surface dynamics in the vicinity of the super-harmonic resonance, we can now experimentally determine their stability regions in the $(\varOmega /\omega _{21}, a_x)$ space. To do so, we fix the forcing amplitude while operating a frequency sweep from high to low frequencies, within the range $\varOmega /\omega _{21} \in [0.45, 0.53]$, by frequency decrements of 10 mHz. Note that a downward frequency sweep ensures recovery of the stability bound between the super-harmonic planar and swirling regimes, as the transition in this direction occurs exactly at the threshold frequency $\varOmega _P(a_x)$ below which the super-harmonic planar motion becomes unstable. On the contrary, since the super-harmonic swirling wave is still stable for frequencies larger than $\varOmega _P(a_x)$ (hysteresis), an upward frequency sweep will maintain the system's response on the swirling branch, thus, it is not suitable to experimentally detect the bifurcation point P. The downward frequency sweep also enables one to detect the bounds that separate the irregular regime from steady planar ($\varOmega = \varOmega _U(a_x)$) and swirling motions ($\varOmega = \varOmega _H(a_x)$).

This procedure is applied for various forcing amplitudes $a_x \in [0.05, 0.4]$, enabling us to build the stability regions diagram displayed on figure 13. All together, the experimental measurements are in very good quantitative agreement with the theoretical regime bounds for $a_x > 0.15$, below which the super-harmonic irregular and swirling regimes appear to be suppressed by dissipative mechanisms, e.g. viscous dissipation occurring in the fluid bulk, sidewall and free-surface boundary layers (Case & Parkinson 1957; Miles 1967; Raynovskyy & Timokha 2020; Bongarzone et al. 2022b) as well as in the neighbourhood of the moving contact line (Keulegan 1959; Dussan 1979; Hocking 1987; Cocciaro, Faetti & Festa 1993; Viola & Gallaire 2018). This last contribution is likely to be important since no particular precautions, such as wall treatment or pre-wetting, have been taken in order to minimize contact angle hysteresis. Below this threshold amplitude, experiments have shown a vanishing super-harmonic contribution to the dynamics, with a harmonic planar motion produced by the SC wave only and well described by the potential linear model, thus suggesting that the DC wave has been entirely killed by the dissipative mechanism. Note that such a suppression of the DC dynamics at low forcing amplitude is reminiscent of what was observed by Reclari et al. (2014) for circular shaking.

Figure 13. Estimates of regime bounds in the $(\varOmega /\omega _{21},a_x)$ plane for a container of diameter $D=0.172\ \text {m}$, filled to a depth $H=1.3$, driven longitudinally and super harmonically at a frequency $\varOmega \approx \omega _{21}/2$: comparison between the theoretical predictions (solid lines) and experimental measurements (markers). Grey thick solid lines: present theoretical predictions. Black empty squares: super-harmonic planar motion. Black crosses: irregular regime. Black filled circles: super-harmonic swirling motion.

6.6. Irregular regime

In this section we provide a more thorough description of the irregular regime. When fixing the forcing frequency slightly below $\omega _{21}/2$ and progressively increasing the forcing amplitude, the free-surface response is first very regular and displays a planar dynamics for low enough forcing amplitudes. Above a threshold amplitude, the dynamics becomes irregular and at large enough amplitudes, the response is again regular, but consists in a swirling motion. Figure 14(a) displays the free-surface response along the vertical middle axis $Y=0$ for increasing forcing amplitudes at a fixed forcing frequency $\varOmega \approx 0.496 \omega _{21}$. The regular regimes (top and bottom panels) are characterized by a constant amplitude of the free-surface oscillations. In contrast, the oscillations of the free surface for intermediate forcing amplitudes (second and third panels) are enclosed into a quasi-periodic envelope, whose frequency linearly increases with the forcing amplitude (see figure 14(b) and Appendix D for the methodology used to compute the envelope). This is very reminiscent of the observations by Royon-Lebeaud et al. (2007) of the irregular regime present in the vicinity of the harmonic resonance under longitudinal forcing. Note that these features are also quantitatively recovered by proceeding with a downwards amplitude sweep. In particular, upward and downward forcing amplitude sweeps provide the same threshold amplitudes between planar and irregular dynamics, and between irregular and swirling regimes. Furthermore, the frequency of the main peak in the power spectrum of the amplitude envelope seems to be a robust feature that does not depend on the sweep direction; see figure 14(b).

Figure 14. (a) Intensity profiles $I(0)$ for various forcing amplitudes and the same forcing frequency $\varOmega \approx 0.496 \omega _{21}$. The oscillations of the free surface are enclosed into an envelope, plotted in black on top of the images. (b) Left: frequency of the main peak in the power spectrum of the envelope (adimensionalized by the forcing frequency $\bar {\varOmega }/2{\rm \pi}$) as a function of the forcing amplitude $a_x$, for the same angular forcing frequency $\varOmega \approx 0.496 \omega _{21}$, and for upwards (white markers) and downwards (black markers) amplitude sweeps. When the envelope is a straight line (as it is the case here for $a_x < 0.15$, which corresponds to a regular planar dynamics of the free surface), the power spectrum is flat and we set the corresponding frequency equal to zero. For $a_x \approx 0.40$, the power spectrum of the envelope is dominated by a low amplitude and small frequency noise, causing a brutal decrease of the ‘burst’ frequency, thus indicating a transition from irregular to regular swirling motion. Right: power spectra of the envelope for $a_x\approx 0.12$, $a_x \approx 0.23$ and $a_x \approx 0.40$. (c) Correlation between $I(R/2)$ and $I(-R/2)$ as a function of time for the same set of forcing parameters as in (a). (d) Left: superposition of $I(R/2)$ (blue) and $I(-R/2)$ (red); right: position of the front contact line $\eta ( R, \pm {\rm \pi}/6, t)$ as a function of time extracted from $I(-R/2)$ (red curve) and from $I(R/2)$ (blue curve). The signals presented in (d) are taken from the full signals used to compute their correlation in (c) for $a_x \approx 0.23$, over the time ranges highlighted in blue and denoted as (1) (maximum of correlation) and (2) (minimum of correlation). (e,f) Normalized power spectra of $\eta (R, -{\rm \pi} /6, t)$ (red curve) and of $\eta (R, {\rm \pi}/6, t)$ (blue curve), and the absolute value of the phase shift between their harmonic and super-harmonic components, where the dynamics of $\eta (R, \pm {\rm \pi}/6, t)$ is considered over (e) time range (1) and (f) time range (2).

To gain more insight on this irregular dynamics, we compute at each time $t_i$ the spatial correlation between $I_{t_i}(R/2)$ and $I_{t_i}(-R/2)$, which we refer to as ${\rm corr} (t_i)$,

(6.1)\begin{equation} {\rm corr}(t_i)=\frac{\displaystyle \sum_n (I_{t_i, n}(R/2)-\bar{I}_{t_i}(R/2))(I_{t_i, n}({-}R/2)-\bar{I}_{t_i}({-}R/2))}{\displaystyle \sqrt{\sum_n (I_{t_i, n}(R/2)-\bar{I}_{t_i}(R/2))^2\sum_n (I_{t_i, n}({-}R/2)-\bar{I}_{t_i}({-}R/2))^2}}, \end{equation}

where $n \in [1, N]$, with $N$ the number of pixels in the vertical direction, and $\bar {I}_{t_i}(y)$ represents the mean of the $N$-element vector $I_{t_i}(y)$. A high and constant correlation is a signature of a steady planar motion, while a low but still constant correlation is characteristic of the steady swirling regime. At intermediary forcing amplitudes – i.e. in the irregular regime – the correlation is a quasi-periodic function of time, with the same quasi-period as the envelope; see figure 14(c).

A comparison between $I(R/2)$ and $I(-R/2)$ on time ranges corresponding to the maximum and minimum of the correlation function reveals that in the time interval where the signals are highly correlated, the motion is planar like (although irregular), while in the time range where they are poorly correlated, the maxima of the right- and left-hand side signals are phase shifted with respect to each other, thus reflecting the presence of a swirling wave; see figure 14(d).

This is further confirmed by the power spectra of the front contact line dynamics along the azimuthal directions $\theta =\pm {\rm \pi}/6$, extracted from $I(R/2)$ and $I(-R/2)$, on time ranges where these signals are highly correlated and where they are poorly correlated; see figure 14(e,f). In both cases, the sloshing wave contains a planar SC wave, as revealed by the vanishing phase shift between the harmonic components of $\eta (R, {\rm \pi}/6, t)$ and $\eta (R, -{\rm \pi} /6, t)$. The wave also contains a super-harmonic component that is responsible for the switching between a planar-like motion (vanishing phase shift between the $\omega _{21}$ components of the $\eta (R, {\rm \pi}/6, t)$ and $\eta (R, -{\rm \pi} /6, t)$ signals, figure 14e) and a swirling dynamics (rotating, symmetry-breaking wave that is super-harmonically oscillating at $\omega \approx \omega _{21}$, figure 14f).

This is again very similar to the features of the irregular regime in the vicinity of the harmonic resonance described by Royon-Lebeaud et al. (2007) that relate the ‘bursts’ in the free-surface oscillation amplitude to the quasi-periodic occurrence of a swirling wave. However, in the case of super-harmonic resonance the irregular regime consists here in the superposition of a stable planar SC wave and of a super-harmonic DC dynamics. The latter is responsible for the irregularity of the total dynamics, by quasi-periodically switching between super-harmonic planar and swirling motions.

6.7. Wave amplitude saturation: theoretical predictions vs experiments

In this last section we provide a more quantitative comparison in terms of wave amplitude saturation between the theoretical predictions according to (5.13) and the experimental measurements. On this point, the dimensional wave amplitude, ${\rm \Delta} \bar {\delta }=\max _{\theta,t}\eta (r=R,\theta,t)- \min _{\theta,t}\eta (r=R,\theta,t)$, is experimentally measured by fixing the forcing amplitude while operating a frequency sweep in two directions. A backward sweep is used so as to follow the right lower planar branch until the sub-critical jump-up transition to swirling ($P$: Poincaré bifurcation) occurs ($\varOmega =\varOmega _{P}(a_x)$). On the other hand, an upward sweep is performed in order to maintain a stable super-harmonic swirling response from bifurcation point $H$ ($\varOmega =\varOmega _{H}(a_x)$) and beyond the threshold frequency $\varOmega =\varOmega _{P}(a_x)$, above which the super-harmonic planar and swirling motions are both stable solutions (right region in the stability chart of figure 13).

For each set of forcing parameters $(a_x, \varOmega /\omega _{21})$, the height in pixel of the wave crest (respectively trough) on the front wall is manually extracted from the corresponding movies and is compared, in the same frame it is extracted from, to the height of the fluid at rest (flagged by a black mark on the container, also used as a scale) to obtain the maximal (respectively minimal) front contact line elevation $\max _{\theta,t}\eta (r=R,\theta,t)$ (respectively $\min _{\theta,t}\eta (r=R,\theta,t)$). This value is converted into metres using the conversion factor provided by the black scale. The resulting amplitude ${\rm \Delta} \bar {\delta }$ is then averaged over three to five cycles of oscillations and finally normalized by the container radius $R$.

The experimental dimensionless wave amplitude ${\rm \Delta} \delta = {\rm \Delta} \bar {\delta }/R$ as a function of the forcing frequency for various forcing amplitudes is displayed in figure 15 together with the theoretical WNL prediction (5.13) (light blue solid lines) and with the linear potential solution (3.4) for comparison (black dashed line).

Figure 15. Quantitative comparison with experimental measurements in terms of finite amplitude saturation for various non-dimensional forcing amplitudes, $a_x$. Black dotted lines: linear potential solution according to (3.4) and (3.5). Light blue solid lines: stable planar and swirling branches predicted by the present WNL model according to (5.13). Markers: experimental measurements. Black empty squares correspond to planar motion whereas black filled circles refer to swirling dynamics. Results are shown for (a) $a_x=0.1628$, (b) $a_x=0.1861$, (c) $a_x=0.2093$, (d) $a_x=0.2325$, (e) $a_x=0.2558$ and (f) $a_x=0.2791$.

The experimental data associated with the two planar branches compare generally well with the present WNL prediction, although the WNL theory slightly underestimates the wave amplitude in the swirling regime. Such an underestimation is not necessarily produced by overlooked nonlinear mechanisms, but it could be instead ascribed to the linear superposition of viscous modes, whose phase shift's characterization would require a dedicated theoretical viscous analysis. The same holds for the asymmetry experimentally observed in the DC wave profile reported in figures 9(c), 11(d) and 12(e).

We recall from § 5 that, at leading order, the wave solution in these two branches is made by the superposition of two planar waves, i.e. an harmonic planar SC component, oscillating in space and time as $\cos {(\varOmega t)}\cos {\theta }$, and a super-harmonic planar DC component, characterized by $\cos {(2\varOmega t+\varPhi )}\cos {2\theta }$, with a phase $\varPhi ={\rm \pi}$ in the left branch and $\varPhi =0$ in the right one. The information on the phase $\varPhi$ is not directly discernible from the amplitude plot of figure 15, but it is contained in the snapshots sequence reported in figure 11(a) for $\varOmega /\omega _{21}<0.5$ and (e) for $\varOmega /\omega _{21}>0.5$. Due to the temporal periodicity of the SC wave, snapshots taken at $t=T/4={\rm \pi} /2\varOmega$ and $t=3T/4=3{\rm \pi} /2\varOmega$ represent temporal nodes for the harmonic component, so that, as a first-order approximation, only the DC component, whose azimuthal spatial structure reads $\cos {({\rm \pi} +\varPhi )}\cos {2\theta }$, is instantaneously left. It is then clear that for $\varOmega /\omega _{21}<0.5$ and $\varPhi ={\rm \pi}$, the free surface maximum is reached at the azimuthal coordinates $\theta =0$ and ${\rm \pi}$, whereas the minimum is at $\theta =\pm {\rm \pi}/2$ (vice versa for $\varOmega /\omega _{21}>0.5$ and $\varPhi =0$). This produces the concave and convex shapes in the instantaneous free surface displayed in figures 11(a) and 11(e), respectively.

Consistently with the stability chart in figure 13 obtained through a backward frequency sweep, the threshold frequency $\varOmega _H(a_x)$, at which the swirling branch becomes stable from lower driving frequency $\varOmega$, is correctly detected. Furthermore, the upward sweep allows us to detect also the jump-down transition from the swirling to the lower right planar branch.

The occurrence of the jump-down transition was to be expected as it is produced by dissipative mechanisms (see also § 6.5), which are overlooked by the present inviscid analysis. The associated damping, which is a function of the wave amplitude and of the forcing acceleration amplitude (see Raynovskyy & Timokha (2018b), Raynovskyy & Timokha (2020) and the discussion in Appendix A of Bongarzone et al. 2022a), is responsible for the modulation in the phase lag between the external driving and the wave response, which was shown by Bäuerlein & Avila (2021) (for unidirectional sloshing waves in a rectangular container) to be of crucial importance for a correct prediction of the jump-down frequency.

The damping coefficient could be tentatively fitted from experiments and phenomenologically introduced a posteriori in amplitude (5.9a) and (5.9b) as done in Appendix A of Bongarzone et al. (2022a). Nevertheless, the jump-down transition in the cases examined in this section (see figure 15) was seen to be extremely sensitive to the frequency sweeping rate. A decrease in the frequency step increment from 5 mHz to 1 mHz (used to produced the swirling branch in figure 15) was observed to give different jump-down frequencies. This is also expected as it is known from the literature that in the multi-solution range, the characteristic of the response mainly depends on the sweep rate (Park, Do & Lopez 2011; Bourquard & Noiray 2019; Yu et al. 2020). Since we did not try frequency increments smaller than 1 mHz, the jump-down frequency predictions as shown in figure 15 are not entirely reliable for fitting the damping at stake in the experiments.

In spite of such limitations, the WNL model is seen to describe fairly well the experimental swirling branch until the measured jump-down frequency. A relatively small departure of the swirling response from the theoretical prediction is typically observed at larger driving amplitude for increasing wave frequency. In agreement with previous studies (Dodge et al. 1965; Ibrahim 2005; Bäuerlein & Avila 2021), our experiments reveal that this is due to the progressive steepening and broadening of the wave crest and troughs, respectively, in the vicinity of the container wall. This nonlinear mechanism eventually becomes strong enough for the WNL model to lose accurateness.

7. Conclusion

In this work the behaviour of sloshing waves in a cylindrical container submitted to longitudinal periodic forcing with driving amplitude $a_x$ and angular frequency $\varOmega$ was investigated. While previous studies of this forcing condition and geometry mostly focused on the investigation of the free-surface response in the vicinity of harmonic resonance, i.e. $\varOmega /\omega _{1n}\approx 1$, the core of the present work was dedicated to the most relevant secondary super-harmonic resonances $\varOmega / \omega _{2n} \approx 1/2$, characterized by the occurrence of a DC dynamics oscillating at a frequency $\omega = 2\varOmega \approx \omega _{2n}$.

Such a super-harmonic resonance was first experimentally observed by Reclari (2013) and Reclari et al. (2014) for circular container motions, but its investigation under different forcing conditions, e.g. longitudinal forcing, seemed to be still unreported.

With the aim to take a further step in this direction, a WNL analysis via the multiple time-scale method together with a dedicated experimental campaign were implemented in order to account for the steady-state free-surface dynamics, and for the symmetry breaking due to the emergence of a DC swirling wave in the vicinity of the super-harmonic resonance.

In similar fashion to Bongarzone et al. (2022a), the WNL analysis was first formalized to tackle the simpler case of harmonic resonances. The outcomes of the model were compared with previous experimental measurements and to former theoretical predictions based on the Narimanov–Moiseev multimodal sloshing theory (Faltinsen et al. 2016; Raynovskyy & Timokha 2020). All together, our analysis addressing the SC wave dynamics was shown to be consistent with the previously reported experimental and theoretical results. In particular, the WNL model successfully captured the regime bounds between SC planar, swirling and irregular waves, and correctly described the close-to-resonance nonlinear behaviour, thus validating the relevance of this theoretical approach.

The WNL analysis was then extended to the more complex case of the super-harmonic resonance. A dedicated lab-scale experiment was set-up to observe and characterize the super-harmonic response to longitudinal forcing. In remarkable agreement with the outcomes of the WNL model, the experimental investigation showed that the free-surface dynamics in the vicinity of the super-harmonic resonance results from the superposition of a permanent, first-order forced harmonic planar SC wave, and of a super-harmonic DC wave that can exhibit either a planar, irregular or swirling dynamics, the latter being responsible for a symmetry breaking in the system's response through equally probable clockwise or anti-clockwise swirling waves. The bounds in the $(a_x, \varOmega /\omega _{21})$ plane between the three different regimes were experimentally retrieved and were shown to be in very good quantitative agreement with the WNL predictions, at least above a threshold forcing amplitude, below which the swirling and irregular dynamics appear to be suppressed by dissipative mechanisms, which are not accounted for by the present inviscid analysis. Finally, the predicted wave amplitude saturation, computed by reconstructing the total flow solution, was compared with the experimentally measured steady-state wave amplitude and was shown to correctly describe the stable planar and swirling branches in the neighbourhood of the super-harmonic resonance.

The fairly good agreement between the theoretical predictions and the experimental findings validates the relevance of the WNL approach to successfully describe the sloshing wave dynamics resulting from nonlinear harmonic and super-harmonic interactions. As discussed in Appendix B, this analysis is not restricted to longitudinal forcing, but can be straightforwardly generalized without any further calculation to any elliptic trajectory, hence recovering the limit of circular sloshing investigated in Bongarzone et al. (2022a). In this respect, the theory of Faltinsen et al. (2016) for elliptical container motions interestingly predicts the occurrence of counter-rotating swirling waves, i.e. propagating in the direction opposed to that of the container motion. The qualitative analogy between the harmonic and super-harmonic system behaviour outlined in this paper would suggest that such counter-propagating swirling waves could also be triggered by exciting the system in the vicinity of the super-harmonic DC resonance, thus calling for new experimental campaigns.