An amplitude equation modelling the double-crest swirling in orbital-shaken cylindrical containers

Abstract

Container motion along a planar circular trajectory at a constant angular velocity, i.e. orbital shaking, is of interest in several industrial applications, e.g. for fermentation processes or in cultivation of stem cells, where good mixing and efficient gas exchange are the main targets. Under these external forcing conditions, the free surface typically exhibits a primary steady-state motion through a single-crest dynamics, whose wave amplitude, as a function of the external forcing parameters, shows a Duffing-like behaviour. However, previous experiments in laboratory-scale cylindrical containers have revealed that, owing to the excitation of super-harmonics, diverse dynamics are observable in certain driving-frequency ranges. Among these super-harmonics, the double-crest dynamics is particularly relevant, as it displays a notably large amplitude response, which is strongly favoured by the spatial structure of the external forcing. In the inviscid limit and with regards to circular cylindrical containers, we formalize here a weakly nonlinear analysis via a multiple-time-scale method of the full hydrodynamic sloshing system, leading to an amplitude equation suitable for describing such a double-crest swirling motion. The weakly nonlinear prediction is shown to be in fairly good agreement with previous experiments described in the literature. Lastly, we discuss how an analogous amplitude equation can be derived by solving asymptotically for the first super-harmonic of the forced Helmholtz–Duffing equation with small nonlinearities.

1. Introduction

Orbital shaking is a method to gently mix the liquid content of a container by its displacement at fixed container orientation along a circular trajectory and at a constant angular velocity. It is used in biological and chemical industrial applications, notably bacterial and cellular cultures (McDaniel & Bailey 1969; Wurm 2004), as an alternative to stirred tanks, where liquid agitation results from a rotating impeller or the rotation of a magnetic rod. In these cultivation protocols, cells are in suspension in an extracellular liquid medium, which serves as buffer for consumables from which they feed and for their secretions. The motion of the liquid prevents sedimentation and homogenizes the concentration of dissolved oxygen and nutrients and of secreted proteins and carbon dioxide. Because of the possible gas exchanges at the free surface, oxygen supply from the container bottom can possibly be circumvented, avoiding the formation of bubbles and thereby the damage that their collapse can exert on cells (Handa-Corrigan, Emery & Spier 1989; Kretzmer & Schügerl 1991; Papoutsakis 1991), sparking interest in the development of large-scale, in the hectolitre range, orbital-shaken bioreactors (Liu & Hong 2001; Jesus et al. 2004; Muller et al. 2007). It is therefore not a surprise that a significant body of research on gas exchange and mixing in these devices has emerged over the last two 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).

Since the shear stresses and, therefore, the mixing are proportional to the velocity gradients in the liquid phase, most of the gas exchange phenomena listed above are directly linked to the liquid motion, with the optimal working conditions essentially dictated by the wave pattern (Reclari 2013). For these reasons, at a more fundamental level, the hydrodynamics of these orbital shaking devices has received recent attention, from both experimental (Reclari et al. 2014; Bouvard, Herreman & Moisy 2017; Moisy, Bouvard & Herreman 2018) and theoretical (Reclari et al. 2014; Horstmann, Herreman & Weier 2020) perspectives, predominantly using linear potential flow models. These models are often complemented with effective viscous damping rates to incorporate the energy dissipation responsible for the phase shifts between wave and shaker, which was also seen to be sometimes responsible for damping-induced symmetry-breaking linear mechanisms resulting in linear spiral wave patterns (Horstmann et al. 2020, 2021). Previous studies mostly made use of classical existing theories for general linear sloshing dynamics, reviewed for instance in Ibrahim (2005) and Faltinsen & Timokha (2009).

In order to refine the linear potential model and, specifically, to predict the occurrence of the super-harmonic wave dynamics observed experimentally (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), Reclari (2013) and Reclari et al. (2014) proposed an inviscid weakly nonlinear (WNL) analysis based on a second-order straightforward asymptotic expansion procedure, which was shown to be capable of capturing the observed resonance frequencies and of characterizing different multiple-crest wave patterns. Among these patterns, the super-harmonic double-crest (DC) wave dynamics is particularly relevant, as it appears to be the most stable and the one that displays the largest amplitude response. However, their analysis, as typical of straightforward asymptotic expansions, suffers from secular terms (Castaing 2005; Nayfeh 2008) and, therefore, it still fails in describing the correct nonlinear behaviour close to both harmonic and super-harmonic resonances.

With regards to the experiments of Reclari (2013) and Reclari et al. (2014), Timokha & Raynovskyy (2017) and Raynovskyy & Timokha (2018a,b) have applied the Narimanov–Moiseev multimodal sloshing theory (Narimanov 1957; Moiseev 1958; Dodge, Kana & Abramson 1965; Faltinsen 1974; Narimanov, Dokuchaev & Lukovsky 1977; Lukovsky 1990). The theory is capable of accurately describing the nonlinear wave dynamics near the primary harmonic resonance, when no secondary resonances occur (Faltinsen, Rognebakke & Timokha 2005; Faltinsen, Lukovsky & Timokha 2016). Despite the fact that the experiments performed by Reclari (2013) and Reclari et al. (2014) were made for non-dimensional fluid depths $H=h/R=1.04$ and $1$, which lie slightly beyond the applicability threshold of the multimodal theory ($H_{th}$ should be $\gtrsim 1.05$ as stated by Raynovskyy & Timokha (2020)) and imposed by the occurrence of secondary resonances, the authors found a quantitative good agreement with the experimental observations associated with hard-spring-type single-crest (SC) swirling.

In the spirit of the aforementioned multimodal theory but with regards to square-base basins, the resonant amplification of higher-order modes for forcing frequency in the vicinity of the primary resonance (secondary or internal resonances) was investigated by Faltinsen et al. (2005), who formalized a so-called adaptive asymptotic modal approach capable of improving the agreement with earlier experiments. A thorough discussion on this regard is also outlined in Chapters 8 and 9 of Faltinsen & Timokha (2009), where the importance of the ratio of tank liquid depth to tank breadth for the occurrence of the internal resonance phenomenon is carefully discussed. Generally speaking, secondary resonance is a broader concept and it may occur even far from the primary resonance zone, as in the case of the DC swirling observed in Reclari et al. (2014). To the knowledge of the authors, the adaptive modal approach has never been extended to super-harmonic system responses of orbital-shaken circular cylindrical containers far from the primary resonance.

For these reasons, it appears that a quantitatively accurate model for the prediction of the diverse wave dynamics observed during the thorough experimental campaign carried out by Reclari (2013) and Reclari et al. (2014) has not been provided yet.

The present work is precisely dedicated to the development of a WNL analysis based on the multiple-time-scale method, which will be seen to successfully capture nonlinear effects for the main additive harmonic resonances as well as the more subtle additive and multiplicative resonance governing the super-harmonic DC swirling. Amplitude equations are rigorously derived in an inviscid framework, which, once amended with an ad hoc damping term as the only tuning parameter, well match the experimental findings of Reclari (2013) and Reclari et al. (2014). Lastly, the obtained amplitude equations for harmonic SC and super-harmonic DC waves are found to be compatible with two well-known one-degree-of-freedom systems: the Duffing and the Helmholtz–Duffing oscillators, respectively.

The paper is organized as follows. The flow configuration and governing equations are introduced in § 2. Section 3 is dedicated to a brief summary of the salient points of the asymptotic model proposed by Reclari et al. (2014), the limitations of which motivated the present work. After tackling the more common case of harmonic SC waves in § 4.1, the WNL amplitude equation governing the super-harmonic DC wave dynamics is derived in § 4.2. Final comments and conclusions are outlined in § 5.

2. Flow configuration and governing equations: potential model

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$. The orbital (circular) shaking motion (see sketch in figure 1) can be represented as the combination of two sinusoidal translations with a ${\rm \pi} /2$ phase shift, thus leading to the following equations of motion for the container axis intersection with the $z=0$ plane, parametrized in cylindrical coordinates ($r$, $\theta$):

(2.1)\begin{equation} \dot{\boldsymbol{X}}_0 = \begin{cases} -\dfrac{d_s}{2}\varOmega_d\sin{\left(\varOmega_d t-\theta\right)}\,\boldsymbol{e}_r,\\ \dfrac{d_s}{2}\varOmega_d\cos{\left(\varOmega_d t-\theta\right)}\,\boldsymbol{e}_{\theta}. \end{cases} \end{equation}

In the classical potential flow limit, i.e. the flow is assumed to be inviscid, irrotational and incompressible, the motion is described in terms of free-surface deformation, $\eta$, and a potential velocity field, $\varPhi _{tot}$, which is typically separated into a container component, $\varPhi _c$, and a fluid component, $\varPhi$. Hence, the liquid motion within the moving container is governed by the Laplace equation,

(2.2)\begin{equation} \Delta\varPhi=\frac{1}{r}\frac{\partial\varPhi}{\partial r}+\frac{\partial^2\varPhi}{\partial r^2}+\frac{1}{r^2}\frac{\partial^2\varPhi}{\partial \theta^2}+\frac{\partial^2\varPhi}{\partial z^2}=0, \end{equation}

subjected to the homogeneous no-penetration condition, $\boldsymbol {\nabla }\varPhi \boldsymbol {\cdot }\boldsymbol {n}=\boldsymbol {0}$, at the solid sidewall and bottom, and by the dynamic and kinematic free-surface boundary conditions at $z=\eta$ (see Ibrahim 2005),

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

(2.3b)\begin{gather} \frac{\partial \eta}{\partial t}+\frac{\partial\varPhi}{\partial r}\frac{\partial \eta}{\partial r}+\frac{1}{r^2}\frac{\partial\varPhi}{\partial \theta}\frac{\partial \eta}{\partial \theta}-\frac{\partial\varPhi}{\partial z}=0, \end{gather}

which have been made non-dimensional by using the container's characteristic length $R$, the characteristic 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. The non-dimensional driving amplitude and angular frequency read $f=d_s\varOmega _d^2/(2g)$ and $\varOmega =\varOmega _d/\sqrt {g/R}$, respectively. When surface tension is accounted for, an additional contact line boundary condition is required at $z=\eta$ and $r=1$, typically written as $\partial \eta /\partial r=\cot {\vartheta }$, where $\vartheta$ is the macroscopic contact angle. Under the classic free-end edge contact line assumption with $\vartheta ={\rm \pi} /2$ adopted here, the latter dynamic equation simply reduces to $\partial \eta /\partial r=0$. This means that the free surface at rest is flat and that a ${\rm \pi} /2$ static contact angle is maintained when the contact line elevation changes dynamically.

Figure 1. Sketch of a cylindrical container of diameter $D=2R$ and filled to a depth $h$. The gravity acceleration is denoted by $g$. The Cartesian inertial reference frame is $O'\boldsymbol {e}'_x\boldsymbol {e}'_y\boldsymbol {e}'_z$, 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. Parameter $d_s$ is the diameter of the circular shaking trajectory, characterized by a driving angular frequency $\varOmega _d$.

3. Linear solution and second-order straightforward asymptotic expansion

In order to enlighten the limitations of the expansion procedure developed by Reclari et al. (2014), which motivates the formalization of the new theoretical framework proposed in the present paper, we briefly recall the salient points. Let us consider the following asymptotic expansion for the flow quantities:

(3.1a)\begin{gather} \varPhi=\varPhi_0+\epsilon\varPhi_1+\epsilon^2\varPhi_2+\text{O}\left(\epsilon^3\right), \end{gather}

(3.1b)\begin{gather} \eta=\eta_0+\epsilon\eta_1+\epsilon^2\eta_2+\text{O}\left(\epsilon^3\right), \end{gather}

together with the further assumption of small driving forcing amplitudes of $\text {O}(\epsilon )$, i.e. $f=\epsilon F$, with $\epsilon$ a small parameter $\epsilon \ll 1$ and the auxiliary variable $F$ of $\text {O}(1)$. Solution $\boldsymbol {q}_0=\{\varPhi _0,\eta _0\}^{\rm T}$ represents the rest state, which has a potential velocity field null everywhere, $\varPhi _0=0$, and a flat interface, $\eta _0=0$, as the contact angle is here assumed to be $\vartheta ={\rm \pi} /2$. Substituting the expansions above in (2.2)–(2.3b), a series of systems at the various orders in $\epsilon$ is obtained. At leading order, (2.2)–(2.3b) reduce to a forced linear system, whose matrix compact form reads

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

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

(3.3a,b)\begin{equation} \boldsymbol{\mathsf{B}} = \begin{pmatrix} 0 & 0 \\ {\rm I}_{\eta} & 0 \end{pmatrix},\quad \boldsymbol{\mathsf{A}} = \begin{pmatrix} \Delta & 0 \\ 0 & -{\rm I}_{\eta}+\dfrac{1}{Bo}\dfrac{\partial\kappa}{\partial\eta} \end{pmatrix}, \end{equation}

where ${\rm c.c.}$ stands for complex conjugate, $\partial \kappa /\partial \eta$ represents the first-order variation of the curvature associated with the small perturbation $\epsilon \eta _1$ and ${\rm 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). In the limit of zero external forcing, i.e. $F=0$, system (3.2) is a linear homogeneous problem which, by seeking for solutions having the following normal form:

(3.4)\begin{equation} \hat{\boldsymbol{q}}_{{mn}}\left(r,z\right)\exp({\text{i}\left(\omega_{mn} t-m\theta\right)})+{\rm c.c.}, \end{equation}

reduces to the classic generalized eigenvalue problem for inviscid capillary–gravity waves:

(3.5)\begin{equation} \left(\text{i}\omega_{mn}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\hat{\boldsymbol{q}}_{{mn}}=\boldsymbol{0}, \end{equation}

where indices $(m,n)$ represent the number of nodal circles and nodal diameters, respectively, with $m$ also commonly known as azimuthal wavenumber. Owing to the normal mode expansion, we note that the operator $\mathcal {A}$ depends on the azimuthal wavenumber, $m$, and, therefore, we denote it by $\mathcal {A}_m$. An exact analytical solution to (3.5) can be readily obtained via a Bessel–Fourier series representation leading to the well-known dispersion relation (Lamb 1993)

(3.6)\begin{equation} \omega_{mn}^2=(k_{mn}+k_{mn}^3/Bo)\tanh{(k_{mn}H)}, \end{equation}

with $H=h/R$ and 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 ${\rm J}'_{m}(k_{mn})=0$.

Despite the existence of this analytical solution, in this work we opt for a numerical scheme based on a discretization technique, where linear operators $\mathcal {B}$ and $\mathcal {A}_m$ are discretized in space by means of a Chebyshev pseudospectral collocation method with a two-dimensional mapping implemented in Matlab, which is analogous to that described by Viola et al. (2018). This numerical technique will enable us to avoid straightforward, but cumbersome calculations, otherwise required in the development of the rest of this work and, particularly, of § 4.2. One must note that when (3.5) is solved numerically as in the present case, additional boundary conditions need to be made explicit in order to regularize the problem on the revolution axis ($r=0$), i.e.

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

(3.7b)\begin{gather} |m|\ge1:\quad \hat{\eta}_{{mn}}=\hat{\varPhi}_{{mn}}=0. \end{gather}

It is also useful to note that owing to the symmetries of the problem, system (3.5) is invariant under the 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}

Convergence of the numerical solution was checked by computing the first 16 modes ($m=0,2,3,4$ with $n=1,2,3,4$), whose corresponding natural frequency values, $\omega _{mn}$, matched the analytical ones given by (3.6) up to the fourth digit for a computational grid $N_r=N_z=60$, with $N_r$ and $N_z$ the number of radial and axial grid points, respectively.

Let us now reintroduce the forcing term on the right-hand side of (3.2). In contradistinction with the cases of unidirectional forcing (Miles 1984a,b), for circular orbits, given the azimuthal periodicity of the associated forcing, the shaking at linear order is expected to excite non-axisymmetric modes only and, specifically, those with $m=1$. Therefore, the linear response to the external forcing can be sought as

(3.9)\begin{equation} \boldsymbol{q}_{1}=F\hat{\boldsymbol{q}}_1^F \exp\left({\text{i}\left(\varOmega t-\theta\right)}\right)+{\rm c.c.}, \end{equation}

with $\hat {\boldsymbol {q}}_1^F$ being the solution of the following forced problem:

(3.10)\begin{equation} \left(\text{i}\varOmega\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{1}\right)\hat{\boldsymbol{q}}_1^F=\boldsymbol{\hat{\mathcal{F}}}_1^F. \end{equation}

The response structure $\hat {\boldsymbol {q}}_1^F$ is here computed numerically, but, in practice, it is formally equivalent to that obtained analytically by Reclari et al. (2014) by projecting the forcing term $\boldsymbol {\hat {\mathcal {F}}}_1$ onto the basis formed by the first-order Bessel functions of the first kind, except that surface tension is retained here because of the finite Bond number. Noting that $\epsilon F=f=d_s\varOmega ^2/(2g)$, in figure 2 the linear solution $\epsilon \boldsymbol {q}_1^F$ from (3.9) is shown (black solid lines) and compared with experimental measurements reported by Reclari et al. (2014) in terms of maximum non-dimensional crest-to-trough contact line amplitudes, $\tilde {\delta }=\delta R/D$, with $\delta (\theta,t)=\eta (r=1,\theta,t)$. Measurements for different values of the non-dimensional shaking diameters, $\tilde {d}_s=d_s/D$, are shown. Blue and green markers in figure 2 correspond to highly nonlinear scenarios manifesting a free-surface breaking, which is therefore subsequently ignored. As discussed by Reclari et al. (2014) and reproduced here, the linear solution describes well the SC wave dynamics for driving frequencies far enough from harmonic resonances and, particularly, for small $\tilde {d}_s$. However, as typical of undamped forced oscillators, the amplitude of the inviscid linear response to the external forcing is proportional to $1/|\varOmega ^2-\omega _{1n}^2|$ and therefore it naturally diverges close to $\omega _{1n}$, thus failing in predicting the close-to-resonance behaviour, e.g. for $\tilde {d}_s=0.02$ at $\varOmega \approx \omega _{11}$. Introduction of viscous dissipation would regularize the divergent behaviour at $\varOmega =\omega _{11}$; however, in the absence of any nonlinear restoring term, the hardening nonlinearity displayed in figure 2 cannot be retrieved.

Figure 2. Markers correspond to the experimentally measured maximum crest-to-trough contact line amplitude (non-dimensional), with $\tilde {\delta }=\delta R/D=\delta /2$, reported by Reclari et al. (2014) for two container diameters, $D=0.144\,\text {m}$ and $D=0.287\,\text {m}$, a non-dimensional depth $\tilde {H}=h/D=0.52$ and five values of $\tilde {d}_s=d_s/D$, as a function of the non-dimensional shaking frequency $\varOmega$ normalized by the natural frequency of the first non-axisymmetric mode, $\omega _{11}=1.3286$ ($m=1$), on the bottom $x$ axis and by that of first non-axisymmetric mode with $m=2$, $\omega _{21}=1.7475$, on the top $x$ axis (the frequency values correspond to $D=0.287\,\text {m}$). Colours denote different wave conditions. Black solid lines: linear potential model solution, from (3.9), computed by solving numerically equation (3.10). Red solid lines: WNL solution close to $\varOmega \approx \omega _{21}/2$, obtained by computing (3.14). Note that in order to ease the comparison with experiments, the non-dimensional $\delta$ was rescaled by a factor $R/D=1/2$, as the container diameter $D$, rather than the container radius $R$, was used by Reclari et al. (2014) to make the equations non-dimensional.

Furthermore, in experiments multiple-crested waves were observed at fractions of the natural frequencies (red markers in figure 2), i.e. the system responses with a frequency which is n times (with n a positive integer) that of the external forcing. Here we refer to such conditions as super-harmonic dynamics (note that the terminology sub-harmonic was used by Reclari et al. (2014) instead). Among these super-harmonics, the DC wave dynamics, occurring at a driving frequency $\varOmega \approx \omega _{21}/2$, was seen to be the most relevant (see figure 2), i.e. the most stable and the one displaying the largest deviation from the linear approximation. This specific multiple-crest dynamics, which is intrinsically nonlinear, is indeed favoured by the azimuthal symmetry of the external forcing. Reclari et al. (2014) tentatively described such a DC dynamics by pursuing the asymptotic analysis up to the second order in $\epsilon$, as in (3.1a) and (3.1b), in order to account for second-order system weak nonlinearities.

At the second order in $\epsilon$, one obtains the following forced linear system:

(3.11)\begin{equation} (\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}})\boldsymbol{q}_2=\boldsymbol{\mathcal{F}}_2=F^2(\boldsymbol{\hat{\mathcal{F}}}_2^{FF} \exp({\text{i}(2\varOmega t-2\theta)})+{\rm c.c.})+F^2\boldsymbol{\hat{\mathcal{F}}}_2^{F\bar{F}}, \end{equation}

where $\boldsymbol {{\mathcal {F}}_2}$ gathers a series of terms produced by the first-order solution through the second-order system nonlinearities. For the sake of brevity, the explicit expression of these forcing terms is here omitted (see Appendix D for details). The overbar denotes the complex conjugate. Also note that amplitude $F$ is actually a real quantity and in the following the superscript ${\bar {F}}$ is used only to indicate forcing terms produced by the combination of the direct and complex conjugate contributions of the first-order response to the external forcing. The right-hand side of (3.11) clearly shows how second-order terms naturally induce a super-harmonic response, whose spatial periodicity is $m=2$, hence precisely corresponding to the DC dynamics experimentally observed. The second forcing term on the right-hand side of (3.11) has $\omega =0$ and $m=0$, i.e. it is steady and axisymmetric. It originates in the leading-order contribution in time- and azimuthal-averaged flow, the so-called mean flow. Equation (3.11) was solved analytically by Reclari et al. (2014) by retaining for convenience only two modes, namely those with $(m,n)=(2,1)$ and $(0,1)$, expected to be the relevant ones. The numerical scheme employed in this work allows us to effortlessly account for all the $(2,n)$ and $(0,n)$ modes simultaneously, as their contribution is directly encompassed in the spatial function $\hat {\boldsymbol {q}}_2^{FF}$ and $\hat {\boldsymbol {q}}_2^{F\bar {F}}$, appearing in the second-order solution,

(3.12)\begin{equation} \boldsymbol{q}_{2}=(F^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(2\varOmega t-2\theta\right)})+{\rm c.c.})+F^2\hat{\boldsymbol{q}}_2^{F\bar{F}}, \end{equation}

whose contributions are computed by solving the following systems:

(3.13)\begin{equation} \left(\text{i}2\varOmega\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2}\right)\hat{\boldsymbol{q}}_2^{FF}=\hat{\boldsymbol{\mathcal{F}}}_2^{FF},\quad -\boldsymbol{\mathsf{A}}_0\hat{\boldsymbol{q}}_2^{F\bar{F}}=\hat{\boldsymbol{\mathcal{F}}}_2^{F\bar{F}}. \end{equation}

The total flow field, obtained through the asymptotic model, is then given by the sum of the first- and second-order solutions:

(3.14)\begin{equation} \boldsymbol{q}=(f\hat{\boldsymbol{q}}_1^F \exp({\text{i}\left(\varOmega t-\theta\right)})+f^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(2\varOmega t-2\theta\right)})+{\rm c.c.})+f^2\hat{\boldsymbol{q}}_2^{F\bar{F}}, \end{equation}

where, in order to eliminate the implicit small parameter $\epsilon$, the amplitudes $\epsilon F$ and $\epsilon ^2 F^2$ are recast in terms of the physical amplitudes $f$ and $f^2$, respectively. The resulting prediction of the maximum crest-to-trough contact line amplitude $\delta (\theta,t)=\eta (r=1,\theta,t)$ is shown in figure 2 for driving frequencies close to $\varOmega /\omega _{21}\approx 0.5$ (see top $x$ axis) as red solid lines. Although this straightforward asymptotic expansion detects the emergence of the super-harmonic DC wave in that frequency window, it completely fails in capturing the correct nonlinear wave amplitude saturation displaying a hardening behaviour clearly visible in figure 2. Once again, the amplitude of the inviscid second-harmonic response is proportional to $1/|\varOmega ^2-\omega _{2n}^2/4|$ and the total solution tends to diverge close to the DC super-harmonic at $\omega _{21}/2$.

Such a symmetric and, in the absence of dissipation, close-to-resonance divergent behaviour is actually expected when performing straightforward asymptotic expansions as they typically suffer from secular (or resonating) terms that must be properly treated (see Castaing (2005) and Nayfeh (2008) among many other references).

4. Weakly nonlinear analysis via multiple-time-scale method

In order to overcome the aforementioned limitations of the straightforward asymptotic expansion procedure and thus to attempt to bridge the gap between theoretical predictions and experimental observations, we conduct in this section a WNL analysis based on the multiple-time-scale method. With the aim of deriving a WNL amplitude equation governing the DC dynamics, we first tackle the simpler problem of SC waves. In both cases we look for a third-order asymptotic solution of the system

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

where the zeroth-order solution, $\boldsymbol {q}_0=\boldsymbol {0}$, is omitted.

4.1. Single-crest dynamics

In § 3 the forcing amplitude $f$ was assumed of order $\epsilon$, thus leading to a linear first-order problem directly forced by external shaking, which produces a divergent response close to harmonic resonances. With regards to SC waves and specifically to the harmonic response at a driving frequency close to that of one of the non-axisymmetric 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. The following analysis is therefore expected to hold for $\varOmega =\omega _{1n}+\lambda$, where $\lambda$ is a small detuning parameter assumed of order $\epsilon ^2$. Lastly, 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 at which the free surface oscillates with angular frequency $\approx \omega _{1n}$. Hence, the following scalings are assumed:

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

with $F$ and $\varLambda$ of $\text {O}(1)$. We note that the forcing amplitude could be assumed of order $\epsilon ^2$ (as the other parameters); however, this complicates unnecessarily the second-order problem without modifying the final amplitude equation, even if the values of its coefficients end up being slightly, up to order $\epsilon$, different.

Although the asymptotic expansion is here pursued up to the third order in $\epsilon$, the procedure of the WNL analysis is essentially equivalent to that of the straightforward asymptotic model discussed in § 3. The major difference lies in the solution form of the leading-order problem that is now a homogeneous problem, as in (3.5). Given the azimuthal periodicity of the external forcing, among all possible natural eigenmodes we assume a leading-order solution as

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

where $\hat {\boldsymbol {q}}_1^{A_1}$ is the eigenmode (computed by solving (3.5)) associated with $(m,n)=(1,n)$ and $\omega _{1n}$ is the corresponding natural frequency (solution of (3.6)).

The complex amplitude $A_1$, a function of the slow time scale $T_2$ and still unknown at this stage of the expansion, describes the slow-time amplitude modulation of the oscillating wave $\hat {\boldsymbol {q}}_1^{A_1}$ and introduces a new arbitrariness in the problem, which must be fixed at a higher order. Eigensurface $\hat {\eta }_1^{A_1}$ and eigenpotential field $\hat {\varPhi }_1^{A_1}$, computed for $\omega _{1n}=\omega _{11}$, are shown in figures 3(a) and 3(b), respectively.

Figure 3. (a–c) Real part of the first-order $\hat {\eta }_1^{A_1}$ and second-order $\hat {\eta }_2^{A_1A_1}$ and $\hat {\eta }_2^{A_1\overline {A_1}}$ free-surface deformations computed for $\omega _{1n}=\omega _{11}$. (d–f) Imaginary part of the associated first-order $\hat {\varPhi }_1^{A_1}$ and second-order $\hat {\varPhi }_2^{A_1A_1}$ and $\hat {\varPhi }_2^{A_1\overline {A_1}}$ potential velocity field. Each response is denoted by its amplitude dependence, $\epsilon A_1$, $\epsilon ^2 A_1\bar {A}_1$ and $\epsilon ^2A_1A_1$. The first-order eigenmode is normalized with the amplitude and phase of the contact line (at $r=1$), such that the free surface $\hat {\eta }_1^{A_1}$ is purely real, whereas $\hat {\varPhi }_1^{A_1}$ is purely imaginary. Note that the second-order mean flow constantly induces an upside-down bell-like axisymmetric interface deformation pushing the free surface downward at the centre of the moving container. Calculations are performed for the case of figure 2, i.e. pure water with $\rho =1000\,\text {m}\,\text {m}^{-3}$, $\gamma =0.072\,\text {N}\,\text {m}^{-1}$, $D=0.287\,\text {m}$ and $\tilde {H}=h/D=0.52$, for which $Bo=2802.8$ and $\omega _{11}=1.3286$.

By pursuing the expansion to the second order, a linear system forced by the first-order solution and analogous to that of (3.11) is obtained (see Reclari (2013) for the full expansion of the original nonlinear governing equation up to the second order). Nevertheless, the forcing terms on the right-hand side are here proportional to $A_1^2$ (super- or second-harmonic) and to $A_1\bar {A}_1$ (mean flow). Thus, we seek for a second-order solution of the form

(4.4)\begin{equation} \boldsymbol{q}_2=A_1\bar{A}_1\hat{\boldsymbol{q}}_2^{A_1\bar{A}_1}+(A_1^2\hat{\boldsymbol{q}}_{2}^{A_1A_1} \exp({\text{i}\left(2\omega_{1n} t-2\theta\right)})+{\rm c.c.}), \end{equation}

with $\hat {\boldsymbol {q}}_2^{A_1\bar {A}_1}$ and $\hat {\boldsymbol {q}}_{2}^{A_1A_1}$ computed numerically and displayed in figures 3(b,d) and 3(c, f), respectively, in terms of second-order free-surface deformations and potential velocity fields evaluated for $\omega _{1n}=\omega _{11}$. From a numerical perspective, we note that the second-order responses can be straightforwardly computed as long as the pairs $(\omega,m)=(2\omega _{1n},2)$ and $(0,0)$ do not correspond to eigenvalues of (3.5), i.e. the second-order operators $(\text {i}2\omega _{1n}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2})$ and $-\boldsymbol{\mathsf{A}}_0$ are non-singular and hence invertible.

With regards to figure 3, it is interesting to note how the second-order mean flow potential velocity field is null everywhere. This can be rigorously proven by first noticing that the mean flow corresponds to a time- and azimuthal-averaged flow, i.e. $\partial /\partial t=\partial /\partial \theta =0$. Moreover, in the inviscid limit, free-surface elevation and potential field have a ${\rm \pi} /2$ phase shift, meaning that the first-order eigenmode can be normalized such that the eigensurface is purely real, whereas the eigenpotential is purely imaginary. Under these conditions, the mean flow forcing term on the right-hand side of the kinematic equation cancels out, so that the associated Laplace equation appears to be constrained by homogeneous Neumann conditions on all the domain boundaries, thus prescribing a trivial constant potential field and therefore a null velocity field. In other words, the second-order mean flow system reduces to forced linear meniscus equation (resulting from (2.3a)) and its conditions at $r=0$ and $r=1$ (both $\partial \hat {\eta }_2^{A_1\bar {A}_1}/\partial r=0$), which prescribes a static mean interface deformation only. Such a result was expected since the second-order mean flow response represents the Eulerian mean flow, which, together with the so-called Stokes drift, contribute to the overall Lagrangian mean flow (see van den Bremer & Breivik (2018) for a thorough review).

While the Stokes drift is a pure kinematic concept, the Eulerian mean flow, often referred to as streaming flow (Bouvard et al. 2017), is generally believed to be of viscous origin, although another appealing interpretation has been recently proposed (Faltinsen & Timokha 2019). Sticking to the well-accepted viscous Eulerian mean flow generation mechanism, it is not a surprise that the absence of viscous boundary layers results in a vanishing Eulerian mean flow.

We now move forward to the $\epsilon ^3$-order problem, which is once again a linear problem forced by combinations of the first- and second-order solutions as well as by the slow time derivative of the leading-order solution and by the external forcing, which was assumed of order $\epsilon ^3$:

(4.5)\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_3=\mathcal{F}_3 ={-}\partial_{T_2} A_1\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_1} \exp({\text{i}\left(\omega_{1n} t-\theta\right)}) + |A_1|^2A_1 \boldsymbol{\hat{\mathcal{F}}}_3^{A_1\bar{A}_1A_1} \exp({\text{i}\left(\omega_{1n} t-\theta\right)})\nonumber\\ &\quad +F\boldsymbol{\hat{\mathcal{F}}}_3^F \exp({\text{i}{\varLambda} T_2})\exp({\text{i}\left(\omega_{1n} t-\theta\right)}) +\text{NRT}+{\rm c.c.}, \end{align}

with $\boldsymbol {\hat {\mathcal {F}}}_3^F=\{0,r/2\}^{\rm T}$ and where $\text {NRT}$ stands for non-resonating terms, which are not relevant for further analysis. As standard in multiple-scale analysis, the indeterminacy introduced by the unknown amplitude $A_1$ is resolved by requiring that secular terms do not appear in the solution to (4.5). Secularity results from all resonant forcing terms in $\mathcal {F}_3$, i.e. all terms sharing the same frequency and wavenumber $(\omega _{1n},1)$ of $\boldsymbol {q}_1$, and in effect all terms explicitly written in (4.5). It follows that a compatibility condition must be enforced through the Fredholm alternative (Friedrichs 2012). Such a compatibility condition imposes the amplitude $B=\epsilon A_1 {\rm e}^{\text {i}{\lambda } t}$ to obey the following normal form:

(4.6)\begin{equation} \frac{{\rm d}B}{{\rm d}t}={-}\text{i}{\lambda} B + \text{i}\,\mu_{{SC}} f + \text{i}\,\nu_{{SC}} |B|^2B, \end{equation}

where the physical time $t=T_2/\epsilon ^2$ has been reintroduced and where forcing amplitude and detuning parameter are recast in terms of their corresponding physical value, $f=\epsilon ^3 F$ and $\lambda =\epsilon ^2\varLambda$. Moreover, the small implicit parameter $\epsilon$ is eliminated by defining the total physical amplitude $A=\epsilon A_1$ (Bongarzone et al. 2021a). The various normal form coefficients, which turn out to be real-valued quantities owing to the absence of dissipation, are computed as scalar products between the adjoint mode, $\hat {\boldsymbol {q}}_1^{A_1 {\dagger} }$, associated with $\hat {\boldsymbol {q}}_1^{A_1}$, and the third-order resonant forcing terms as follows:

(4.7a)\begin{gather} \text{i}\,\mu_{{SC}}=\frac{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\boldsymbol{\mathsf{B}}\hat{\mathcal{F}}_3^F\rangle}{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\mathcal{B}\hat{\boldsymbol{q}}_1^{A_1}\rangle}=\frac{\int_{z=0}^{}r\bar{\hat{\eta}}_1^{A_1 {\dagger}}/2\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\varPhi}_1^{A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\eta}_1^{A_1})\, r\text{d}r}, \end{gather}

(4.7b)\begin{gather} \text{i}\,\nu_{{SC}}=\frac{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\boldsymbol{\mathsf{B}}\boldsymbol{\hat{\mathcal{F}}}_3^{A_1\bar{A}_1A_1}\rangle}{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\mathcal{B}\hat{\boldsymbol{q}}_1^{A_1}\rangle}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_1\bar{A}_1A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_1\bar{A}_1A_1})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\varPhi}_1^{A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\eta}_1^{A_1})\, r\text{d}r}.\end{gather}

Here $\hat {\boldsymbol {q}}_1^{A_1 {\dagger} }=\bar {\hat {\boldsymbol {q}}}_1^{A_1}$, since the inviscid problem is self-adjoint with respect to the Hermitian scalar product $\langle \boldsymbol {a},\boldsymbol {b}\rangle =\int _{\varSigma }^{}\bar {\boldsymbol {a}}\boldsymbol {\cdot }\boldsymbol {b}\,\text {d}V$, with $\boldsymbol {a}$ and $\boldsymbol {b}$ two generic vectors (see Viola et al. (2018) for a thorough discussion and derivation of the adjoint problem). For the sake of brevity, the explicit expression of $\boldsymbol {\hat {\mathcal {F}}}_3^{A_1\bar {A}_1A_1}$ is omitted, as it only involves straightforward calculations, i.e. a Taylor expansion of nonlinear governing equations and boundary conditions (2.2)–(2.3b) around the rest state $\boldsymbol {q}_0=\boldsymbol {0}$. Here we simply denote with the subscripts ‘${dyn}$’ and ‘${kin}$’ the forcing components appearing in the dynamic and kinematic boundary conditions, respectively.

By turning to polar coordinates, $B=|B|{\rm e}^{\text {i}\varTheta }$, splitting the modulus and phase parts of (4.6) and looking for stationary solution, ${\rm d}/{\rm d}t=0$ with $|B|\ne 0$, the following implicit relation is obtained:

(4.8)\begin{equation} \tilde{d}_s\varOmega^2\mp \frac{\left({\lambda}-\nu_{{SC}}|B|^2\right)|B|}{\mu_{{SC}}}=0, \end{equation}

or, in a more common polynomial form,

(4.9)\begin{equation} P\left(|B|\right)=|B|^3-\frac{{\lambda}}{\nu_{{SC}}}|B|\pm\frac{\mu_{{SC}}\tilde{d}_s\varOmega^2}{\nu_{{SC}}}=0, \end{equation}

where $f=\tilde {d}_s\varOmega ^2$, $\lambda =(\varOmega -\omega _{1n})$ and the $\mp$ signs correspond to the phases $\varTheta =0$ and ${\rm \pi}$, respectively. The two branches prescribed by (4.8) for $|B|$ as a function of $\varOmega$ at a fixed non-dimensional shaking diameter $\tilde {d}_s$ can be easily computed using the Matlab function fimplicit. After evaluating the stable and unstable stationary solutions for $|B|$ and $\varTheta$, the total SC wave solution is reconstructed as

(4.10)\begin{equation} \boldsymbol{q}_{SC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2. \end{equation}

4.1.1. Experiments versus WNL prediction: wave amplitude

In figure 4 the WNL prediction in terms of maximum crest-to-trough contact line amplitude, $\Delta \tilde {\delta }$, for SC waves is compared with two sets of experimental measurements and with the potential linear solution (3.9). In comparison with the linear theory presented in § 3, the agreement with experiments improves for different shaking diameters and for different harmonic resonances, e.g. those associated with modes $(m,n)=(1,1)$ and $(1,2)$ of figure 4(a). The hardening nonlinearity is correctly captured and the amplitude prediction matches well the measurements until the free surface eventually breaks and the wave regime leaves the WNL regime, hence suggesting the little relevance of dissipative effects attributable to viscosity in this regime.

Figure 4. Comparison of the WNL prediction for SC waves with experiments in terms of maximum crest-to-trough contact line amplitude (non-dimensional), $\Delta \tilde {\delta }$. (a) Black pentagons correspond to the experimental measurements presented in figure 4.30 of Reclari (2013) (R13) for $\tilde {d}_s=d_s/D=0.01$, where the first two non-axisymmetric modes $(m,n)=(1,1)$ and $(1,2)$ are detected. Dotted black lines: solution of the linear potential model according to (3.9). Light blue and blue lines: WNL SC prediction (4.10). Unstable branches are represented as a dashed line. Violet solid lines: theoretical prediction by Raynovskyy & Timokha (2018a) (R&T18) (see § 4.1.2 for further comments). (b) Same as (a) with the black filled circles corresponding to the measurements of Reclari et al. (2014) (R14) reported in figure 2 for $\tilde {d}_s=0.02$. The values of normal form coefficients $\mu _{SC}$ and $\nu _{SC}$ computed for $\varOmega \approx \omega _{11}$ and $\omega _{12}$ (see bottom and top $x$ axes) are given in table 1, together with the corresponding values of $\tilde {H}$, $D$ and natural frequencies $\omega _{mn}$ used in this calculation.

Table 1. Values of the amplitude equation coefficients $\mu _{SC}$ and $\nu _{SC}$ used to produce figure 4.

However, one must note that in this WNL approach the driving frequency is essentially fixed around that of a unique non-axisymmetric natural mode, $\varOmega \approx \omega _{1n}$. Consequently, when performing the analysis for a mode $(1,n)$, the influence of all other modes is completely overlooked. In consequence, the accuracy of the asymptotic solution rapidly deteriorates moving away from harmonic resonances, when compared with the linear solution (3.9), which turns out to be more accurate. This is visible looking at the bottom stable branch in the multi-solution range of figure 4(b) or by looking at the driving frequency window $\varOmega \in [0.7\omega _{12},0.9\omega _{12}]$ in figure 4(a). In other words, the detuning parameter should be small in order for the present WNL analysis, based on a single-mode expansion, to hold. On this regard, as no other natural frequencies are encountered for $\varOmega <\omega _{11}$, an exception is made for the left-hand branch associated with the harmonic resonance of the first (or fundamental) non-axisymmetric mode, where an excellent agreement of the single-mode prediction, comparable to that of the linear solution, lasts until $\varOmega \approx 0$.

4.1.2. Comparison with the multimodal theory of Raynovskyy & Timokha (2018a)

The violet solid curves reported in figure 4 correspond to the predictions associated with the $\omega _{11}$ SC swirling from the Narimanov–Moiseev multimodal sloshing theory employed by Raynovskyy & Timokha (2018a) (R&T18) (only the stable branches are reported). Their curves have been here carefully reproduced by manually sampling those reported in figure 3 of R&T18 in the range of frequency available, i.e. $\varOmega /\omega _{11}\in [0.8,1.3]$. By looking at the upper stable branch, although an increasing discrepancy is observed for increasing wave amplitudes, one can see that both analyses are in fairly good agreement with experiments and with each other until wave breaking eventually occurs. Such a discrepancy could be tentatively attributed to the different definition of the detuning parameter employed in R&T18, i.e. $\epsilon ^2\varLambda _{R\&T18}=\omega _{1n}^2/\varOmega ^2-1$. On the other hand, by looking at the lower stable branch, one sees that, whereas the jump-up frequencies according to R&T18 and to the present model essentially coincide, the discrepancy between the two predictions increases at larger frequency, i.e. $\varOmega /\omega _{11}>1$, with that of R&T18 being closer to the linear potential model. One should also note that, in contradistinction with our analysis, that of R&T18 accounts for damping and predicts the jump-down transition visible in figure 4(a). This damping was essentially fitted from the experimental measurements and, specifically, from the jump-down frequency occurring at a larger frequency, once the wave-breaking regime is fully established, i.e. $\varOmega /\omega _{11}=1.27$ for $\tilde {d}_s=0.01$ and $\varOmega /\omega _{11}=1.45$ for $\tilde {d}_s=0.02$ (see figure 4). However, experiments suggest that the damping effect on the curves displayed in figure 4 would not be easy to observe, even for $\tilde {d}_s=0.01$. Indeed, the motion undergoes a SC wave breaking, thus entering in a fully nonlinear regime, where both our analysis and that of R&T18 loose predictive power. We therefore decided to discard damping while comparing our results with the close-to-harmonic resonance experiments from Reclari (2013) and Reclari et al. (2014).

4.1.3. The Duffing oscillator analogy

Mass–spring models are widely employed in several engineering fields, e.g. in aerospace engineering, for the description of close-to-resonance sloshing motions (Moiseev 1958; Bauer 1966; Dodge 2000), where nonlinearities are of crucial importance. The most popular driven nonlinear mass–spring model is that developed by Duffing (1918), who added a cubic nonlinear spring deformation (cubic term) to the classical driven harmonic oscillator:

(4.11)\begin{equation} \ddot{x}+2\sigma\dot{x}+ x+c_3 x^3 = p\cos{\varOmega t}, \end{equation}

where $\sigma$ is the damping coefficient and where, depending on the sign of $c_3$, the resonance curve bends and the nonlinear resonance frequency shifts, i.e. it decreases for a softening spring ($c_3 < 0$), whereas it increases for a hardening spring ($c_3 > 0$), thus explaining the original observation of Duffing on vibration mechanism. Ockendon & Ockendon (1973) showed via asymptotic expansion of the potential flow solution in the neighbourhood of a harmonic resonance that for small external forcing amplitudes, sloshing in a two-dimensional rectangular container responds exactly as an undamped Duffing oscillator (with $\sigma =0$). In Appendix B, we briefly show that, as expected, the same holds for close-to-harmonic-resonance sloshing in orbital-shaken cylindrical containers, whose formal amplitude equation, starting from the full inviscid hydrodynamic system, was derived in § 4.1 (see (4.6)). Typically, when the Duffing equation is employed to model close-to-resonance responses in sloshing dynamics and experimental measurements are available, the nonlinear coefficient is often computed by fitting the experimental measurements. Recently, with regards to quasi-two-dimensional rectangular containers laterally excited, Bäuerlein & Avila (2021) have carried out careful quantitative comparisons between experiments and theoretical predictions from the damped Duffing equation, showing that their actual sloshing system is remarkably well described by the forced-damped Duffing oscillator. Nevertheless, for increasing wave amplitude responses, experiments deviate from the Duffing solution, which is not capable of predicting correctly the phase lag between driving and response, shown to be the key factor for an accurate estimation of the sloshing amplitude of the maximal nonlinear resonance (Cenedese & Haller 2020; Bäuerlein & Avila 2021). We note that, by analogy with the undamped Duffing equation, the WNL analysis formalized in § 4.1 exacerbates this aspect, since, owing to the lack of dissipation, it can only predict the classic phase-lag bounds, $0$ and ${\rm \pi}$ (see Appendix A for further comments in this regard). Nevertheless, one should notice that this intrinsic limitation turns out to be unimportant in cases as those of figure 4, where for increasing amplitude a wave breaking eventually occurs and the WNL theory as well as the Duffing mechanical analogy no longer apply.

4.2. Double-crest dynamics

We now tackle the DC wave response. Its formalization is slightly more subtle, as it requires a new reordering of the small control parameter magnitudes as well as an unusual form of the leading-order problem, involving both a homogeneous and a particular solution. We recall that the DC dynamics in figure 2 occurs at a driving frequency $\varOmega \approx \omega _{21}/2$. Results at the end of this section are presented for mode $(2,1)$, for which experiments are available; however, for the sake of generality, we formalize the analysis for any mode $(2,n)$, i.e. $\varOmega =\omega _{2n}/2+\lambda$, where $\lambda$ is the small detuning parameter.

4.2.1. Formalism

To determine a suitable scaling for the forcing amplitude $f$ and detuning parameter $\lambda$ it is instructive to look at the experimental measurements shown in figure 2 for $\varOmega$ close to $\omega _{21}/2$. One can see that approaching $\varOmega \approx \omega _{21}/2$ from lower frequencies, the DC wave emerges on the top of a SC dynamics, with the latter being correctly described by the linear solution, which still behaves well as $\omega _{21}/2$ is far enough from the primary harmonic resonance occurring at $\omega _{11}$. It follows that the forcing amplitude $f$ and detuning $\lambda$ could be retained at order $\epsilon$ and the first-order problem takes the form (3.2), with $\boldsymbol {\mathcal {F}}_1=F\{0,r/2\}^{\rm T}{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}=F\boldsymbol {\hat {\mathcal {F}}}_1^F{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}$, with $f=\epsilon F$ and $\lambda =\epsilon \varLambda$.

Furthermore, in § 3 we have shown how, close enough to the super-harmonic resonance, the divergent behaviour is produced by a second-order resonating term, which breaks the straightforward expansion, as $\epsilon ^2$-order terms should not become larger than the $\epsilon$-order ones. In the following, this asymptotic breakdown is overcome by assuming that the leading-order solution is given by the sum of (i) a particular solution, given by the linear response to the external forcing computed by solving (3.10), and (ii) a homogeneous solution, represented by the natural mode $(m,n)=(2,n)$, obtained by solving the generalized eigenvalue problem (3.5), up to an amplitude to be determined at higher orders. The second-order resonating term will then require, in the spirit of multiple-time-scale analysis, an additional second-order solvability condition, complementing the third-order non-resonance condition already obtained in the SC wave WNL model. This suggests that two slow time scales exist, namely $T_1$ and $T_2$, with $T_1$ one $\epsilon$-order faster than $T_2$, hence implying that quadratic nonlinearities are stronger than cubic ones. To summarize, the fundamental scalings underpinning the WNL multiple-scale expansion for DC waves are the following:

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

The first-order solution reads

(4.13)\begin{equation} \boldsymbol{q}_{1}=A_2\left(T_1,T_2\right)\hat{\boldsymbol{q}}_1^{A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})+F\hat{\boldsymbol{q}}_1^F \exp({\text{i}\left(\left(\omega_{2n}/2\right) t-\theta\right)}) {{\rm e}^{\text{i}\varLambda T_1}}+{\rm c.c.}, \end{equation}

where the unknown slow-time amplitude modulation, $A_2$, is here a function of the two time scales $T_1$ and $T_2$, while the amplitude of the particular solution simply equals the forcing amplitude and $\hat {\boldsymbol {q}}_1^F$ is computed from (3.10) for $\varOmega =\omega _{2n}/2$.

The second-order linearized forced problem reads

(4.14)\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_2={\boldsymbol{\mathcal{F}}_2}={\boldsymbol{\mathcal{F}}_2}^{ij}-\frac{\partial A_2}{\partial T_1}\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2}\exp({\text{i}\left(\omega_{2n} t-2\theta\right)}) \nonumber\\ &\qquad - \text{i}{\varLambda} F\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{{F}} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}) {\exp({\text{i}\varLambda T_1})}+{\rm c.c.} \end{align}

The first-order solution is made up of four different contributions of amplitude $A_2$, $\bar {A}_2$, $F$ and $\bar {F}$; therefore it generates 10 different second-order forcing terms, here denoted by $\boldsymbol {\mathcal {F}}_2^{ij}$, which exhibit a certain frequency and azimuthal periodicity $(\breve {\omega },\breve {m})$.

The additional two forcing terms stem from the time derivative of the first-order solution (4.13) with respect to the first-order slow time scale $T_1$. In order to interpret the last term in (4.14), it is worth first noting that, while the amplitude of the linear solution (3.9), computed at a generic driving frequency, grows with $\varOmega$ as $F/|\varOmega ^2-\omega _{11}^2|=\tilde {d}_s\varOmega ^2/|\varOmega ^2-\omega _{11}^2|\propto \varOmega ^2/|\varOmega ^2-\omega _{11}^2|$, in the WNL model for DC waves, the amplitude of the particular solution (4.13) is proportional to $F/|\omega _{21}^2/4-\omega _{11}^2|=\tilde {d}_s\varOmega ^2/|\omega _{21}^2/4-\omega _{11}^2|\propto \varOmega ^2$, since the driving frequency was frozen at $\varOmega =\omega _{21}/2+\lambda$, with the small detuning parameter $\lambda$ contributing to modify its phase, but not its amplitude. This leads to an increasing discrepancy between (3.9) and the leading-order particular solution (4.13) away from the super-harmonic resonance. The response to the forcing term proportional to $\varLambda F$ in (4.14) can be then interpreted as a second-order correction of the amplitude of the first-order particular solution accounting for a detuning from the exact resonance through $\varLambda F\propto \tilde {d}_s\varOmega ^2(\varOmega -\omega _{2n}/2)$ and contributing to improve the asymptotic approximation in a wider range of driving frequency in the neighbourhood of the super-harmonic frequency.

None of the forcing terms in (4.14) is resonant, as their oscillation frequency and azimuthal wavenumber differ from those of the leading-order homogeneous solution, except the term produced by the second harmonic of the leading-order particular solution, i.e. $\boldsymbol {\mathcal {F}}_2^{FF}=F^2\boldsymbol {\hat {\mathcal {F}}}_2^{FF} \exp ({\text {i}(\omega _{2n}-2\theta )}){\exp ({\text {i}2\varLambda T_1})}+{\rm c.c.}$ To avoid secular terms, a second-order compatibility condition is imposed, requiring that the following normal form is verified:

(4.15)\begin{equation} \frac{\partial A_2}{\partial T_1}=\text{i}\,\mu_{{DC}} F^2{\exp({\text{i}2\varLambda T_1})}, \end{equation}

with $\mu _{{DC}}$ computed as before, i.e.

(4.16)\begin{equation} \text{i}\,\mu_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{2_{{dyn}}}^{FF}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{2_{{kin}}}^{FF})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}. \end{equation}

Taken alone, the dynamics resulting from (4.15) is, however, of little relevance, since the solution, i.e. the frequency–response curve, would still diverge (symmetrically) to infinity for $\Delta =\varOmega -\omega _{2n}/2\rightarrow 0$ 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$. This means that the expansion must be pursued up to the following order, and thereby that we must solve for the second-order solution.

By substituting (4.15) in the forcing expression, (4.14) can be rearranged as follows:

(4.17)\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_2= \boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+\boldsymbol{\mathcal{F}}_{2_{RT}}^{i,j}\nonumber\\ &\quad =\boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+F^2(\boldsymbol{\hat{\mathcal{F}}}_2^{FF}-\text{i}\,\mu_{{DC}} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2}) \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})\,{\exp({\text{i}2\varLambda T_1})}+{\rm c.c.}, \end{align}

where the subscripts ‘${NRT}$’ and ‘${RT}$’ denote non-resonating (whose frequencies and azimuthal periodicities are gathered in table 2) and resonating terms, respectively. Note that the term proportional to ${\varLambda } F$ has been included in the non-resonating forcing terms, whereas the resonant term is written explicitly. The compatibility condition is now trivially satisfied, meaning that the new forcing term is orthogonal to the adjoint mode, $\hat {\boldsymbol {q}}_1^{A_2{\dagger} }=\bar {\hat {\boldsymbol {q}}}_1^{A_2}$, by construction and therefore, according to the Fredholm alternative, a non-trivial solution exists. Hence, we seek for a second-order solution having the following form:

(4.18)\begin{align} \boldsymbol{q}_2&=A_2\bar{A}_2\hat{\boldsymbol{q}}_2^{A_2{\bar{A}_2}}+F^2\hat{\boldsymbol{q}}_2^{F\bar{F}}\nonumber\\ &\quad +A_2^2\hat{\boldsymbol{q}}_2^{A_2A_2} \exp({\text{i}\left(2\omega_{2n}t-4\theta\right)})+{\varLambda} F\hat{\boldsymbol{q}}_2^{{\varLambda} F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}){\exp({\text{i}\varLambda T_1})}\nonumber\\ &\quad +A_2F\hat{\boldsymbol{q}}_2^{A_2F}\exp({\text{i}\left(\left(3\omega_{2n}/2\right)t-3\theta\right)}){\exp({\text{i}\varLambda T_1})}\nonumber\\ &\quad +A_2\bar{F}\hat{\boldsymbol{q}}_2^{A_2\bar{F}} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}){\exp({-\text{i}\varLambda T_1})}\nonumber\\ &\quad +F^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(\omega_{2n}t-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}+{\rm c.c.} \end{align}

All non-resonant responses in (4.18) are handled similarly, i.e. they are computed in Matlab by performing a simple matrix inversion using standard LU solvers (as in §§ 3 and 4.1). As anticipated above, 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 (4.16) ensures that a non-trivial solution for $\hat {\boldsymbol {q}}_2^{FF}$ exists. Diverse approaches can be followed to compute this response. Here such a response is computed by using the pseudoinverse matrix of the singular operator (Orchini, Rigas & Juniper 2016). Another possible approach is given in Appendix A of Meliga, Gallaire & Chomaz (2012), where a two-step regularization procedure, involving an intermediate factious solution for $\hat {\boldsymbol {q}}_2^{FF}$, is employed. We also note that in (4.18), exactly as in (4.4), a second-order homogeneous solution has not been accounted for as its introduction would be irrelevant to the final solution.

Table 2. First-order linear solutions and second-order non-resonating forcing terms gathered by their amplitude dependency and corresponding azimuthal and temporal periodicity $(\breve {m},\breve {\omega })$. Six terms have been omitted as they are the complex conjugates of $\epsilon A_2$, $\epsilon F$, $\epsilon ^2 A_2A_2$, $\epsilon ^2\varLambda F$, $\epsilon ^2A_2F$ and $\epsilon ^2 A_2\bar {F}$.

The first-order solutions together with all the non-resonating second-order responses are shown in the various panels of figure 5, where the two leading-order contributions, $\epsilon A_2$ and $\epsilon F$, corresponding to the DC and SC waves, respectively, can be identified. Moreover, we note that the second-order response proportional to $\epsilon ^2{\varLambda } F$ has a spatial structure similar to that of the leading-order response $\epsilon F$, as it represents the second-order correction to the latter caused by small frequency shifts of order $\epsilon$.

Figure 5. Real part of the first-order (a) $\hat {\eta }_1^{A_2}$ and (e) $\hat {\eta }_1^{F}$ and non-resonating second-order (b) $\hat {\eta }_2^{A_1\bar {A}_2}$, (c) $\hat {\eta }_2^{A_2A_2}$, (d) $\hat {\xi }_2^{\varLambda _1 F}$, ( f) $\hat {\eta }_2^{F\bar {F}}$, (g) $\hat {\eta }_2^{A_2F}$ and (h) $\hat {\eta }_2^{A_2\bar {F}}$ free-surface deformations computed for $\omega _{2n}=\omega _{21}$. The first-order eigenmode is normalized with the amplitude and phase of the contact line (at $r=1$), such that the free surface $\hat {\eta }_1^{A_2}$ is purely real, whereas $\hat {\varPhi }_1^{A_2}$ is purely imaginary. Note that the second-order mean flow $\hat {\eta }_2^{F\bar {F}}$ constantly induces an upside-down bell-like axisymmetric interface deformation pushing the free surface downward at the centre of the moving container, by analogy with the effect produced by $\hat {\eta }_2^{A_1\bar {A}_1}$ for SC waves, as the two responses are essentially equivalent up to a pre-factor. Here the mean flow $\hat {\eta }_2^{A_2\overline {A_2}}$ for DC waves pushes the interface upward at the wall (same as $\hat {\eta }_2^{F\bar {F}}$) and, at the same time, downward in an annular region at intermediate radial coordinates, without altering the free-surface elevation at the container revolutions axis.

Lastly, at third order in $\epsilon$, the problem reads

(4.19)\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_3=\boldsymbol{\mathcal{F}}_3\nonumber\\ &\quad ={-}\frac{\partial A_2}{\partial T_2 }\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)}) - \text{i}2{\varLambda} F^2\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(\omega_{2n}-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}\nonumber\\ &\qquad + |A_2|^2A_2 \boldsymbol{\hat{\mathcal{F}}}_3^{A_2\bar{A}_2A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})+A_2F^2\boldsymbol{\hat{\mathcal{F}}}_3^{A_2F\bar{F}} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})\nonumber\\ &\qquad +{\varLambda} F^2\boldsymbol{\hat{\mathcal{F}}}_3^{\varLambda FF} \exp({\text{i}\left(\omega_{2n}{t}-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}+\text{NRT}+{\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 for the full expression of $\boldsymbol {\mathcal {F}}_2$ and $\boldsymbol {\mathcal {F}}_3$). By noticing that the second and last forcing terms share the same amplitude dependence, i.e. ${\varLambda } F^2$, they can be recast into a single forcing term, say ${\varLambda } F^2\hat {\boldsymbol {\mathcal {F}}}_3^{{\varLambda } FF} \exp ({\text {i}(\omega _{2n}{t}-2\theta )}){\exp ({\text {i}2\varLambda T_1})}+{\rm c.c.}$

Once again, all terms explicitly written in (4.19) are resonant, as they share the same pair $(\omega _{2n},2)$ as the first-order homogeneous solution; hence a third-order compatibility condition, leading to the following normal form, must be enforced:

(4.20)\begin{equation} \frac{\partial A_2}{\partial T_2}=\text{i}\zeta_{{DC}} {\varLambda} F^2{ \exp({\text{i}2\varLambda T_1})}+\text{i}\,\chi_{{DC}}A_2F^2+\text{i}\,\nu_{{DC}} |A_2|^2A_2, \end{equation}

with

(4.21a)\begin{gather} \text{i}\,\zeta_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{{\varLambda} FF}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{{\varLambda} FF})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}, \end{gather}

(4.21b)\begin{gather} \text{i}\,\chi_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_2F\bar{F}}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_2F\bar{F}})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}, \end{gather}

(4.21c)\begin{gather} \text{i}\,\nu_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_2\bar{A}_2A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_2\bar{A}_2A_2})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}. \end{gather}

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 (4.15) and (4.20) into a single equation 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 amplitude $A=\epsilon A_2$. This is achieved by summing (4.15) and (4.20) along with their respective weights $\epsilon ^2$ and $\epsilon ^3$, thus obtaining

(4.22)\begin{equation} \frac{{\rm d}B}{{\rm d}t}={-}\text{i}\,(2{\lambda}-\chi_{{DC}} f^2) B + \text{i}\,(\zeta_{{DC}}{\lambda} + \mu_{{DC}}) f^2 + \text{i}\,\nu_{{DC}} |B|^2B, \end{equation}

where the change of variable $A=B{\rm e}^{\text {i}2 {\lambda } t}$ has been introduced for convenience. As in § 4.1, by turning to polar coordinates, $B=|B| {\rm e}^{\text {i}\varTheta }$, splitting the modulus and phase parts of (4.22) and looking for a stationary solution, ${\rm d}/{\rm d}t=0$ with $|B|\ne 0$, the following implicit relation is obtained:

(4.23)\begin{equation} \tilde{d}_s\varOmega^2-\sqrt{\left(2{\lambda}-\nu_{{DC}}|B|^2\right)|B|/\left[\chi_{{DC}}|B|\pm \left(\zeta_{{DC}}{\lambda}+\mu_{{DC}}\right)\right]}=0, \end{equation}

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

Although two more terms appear in (4.22) and the dependence on the forcing amplitude is different with respect to the SC case, i.e. $f^2$ instead of $f$, thus leading to the square root in (4.23), amplitude equation (4.22) is reminiscent of that given in (4.6). Indeed, (4.22) contains essentially three contributions,

(4.24a–c)\begin{equation} \lambda \leftrightarrow \left(2{\lambda}-\chi_{{DC}}f^2\right),\quad \mu_{{SC}}f\leftrightarrow\left(\zeta_{{DC}}{\lambda}+\mu_{{DC}}\right)f^2,\quad \nu_{{SC}}\leftrightarrow\nu_{{DC}}, \end{equation}

in order: a detuning term (forcing-amplitude-dependent), an additive (quadratic) forcing term (forcing-frequency-dependent) and the classic cubic restoring term. Hence, the same qualitative hardening or softening nonlinear behaviours as well as hysteresis, typical features of the Duffing equation, are expected under the hypotheses of the present analysis.

The total flow solution predicted by the WNL model for DC waves and reconstructed as

(4.25)\begin{equation} \boldsymbol{q}_{DC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2 \end{equation}

is compared in figures 6 and 7 with experiments from Reclari (2013) and Reclari et al. (2014) (see also figure 2).

Figure 6. (a–d) The WNL prediction for DC waves versus experiments by Reclari et al. (2014) (R14) (reported in figure 2) in terms of $\Delta \tilde {\delta }$ for a container diameter $D=0.144\,\text {m}$, $\tilde {H}=0.52$ and for various $\tilde {d}_s$. Dotted black lines: linear model (3.9). Red lines: straightforward asymptotic (3.14). Light-blue solid and dashed lines: stable and unstable branches, respectively, predicted by the WNL solution (4.25). The normal form coefficient values for this configurations are $\chi _{{DC}}=2.755$, $\zeta _{{DC}}=0.150$, $\mu _{{DC}}=0.129$ and $\nu _{{DC}}=10.018$. Violet lines: WNL solution (4.28) including the $\epsilon ^3$-order correction discussed in § 4.2.3.

Figure 7. Same as figure 6 (same colour convention) but for a container diameter $D=0.287\,\text {m}$ and $\tilde {H}=0.50$ (experiments are reported in figure 4.19 of Reclari (2013) (R13)). Normal form coefficients: $\chi _{{DC}}=2.709$, $\zeta _{{DC}}=0.169$, $\mu _{{DC}}=0.134$ and $\nu _{{DC}}=9.885$.

4.2.2. Experiments versus WNL prediction: wave amplitude

In figures 6 and 7, the WNL prediction of DC waves according to (4.25) (light-blue solid and dashed lines) is quantitatively compared with the experimental measurements from Reclari (2013) and Reclari et al. (2014) in terms of maximum non-dimensional crest-to-trough contact line amplitude $\Delta \tilde {\delta }$ for different values of the shaking diameters $\tilde {d}_s$ corresponding to those of figure 2 in the frequency window close to $\omega _{21}/2$.

The improvement gained through the formal WNL analysis, when compared with the linear and straightforward asymptotic models, is striking. The amplitude equation model correctly predicts the surge of the DC swirling and the resulting finite-amplitude saturation via hardening nonlinear mechanism, thus markedly narrowing the gap with experiments for all values of $\tilde {d}_s$ considered and for different container configurations.

Notwithstanding such an improvement, figure 6 highlights the main limitation of the present amplitude equation model for DC waves. Indeed, one notices that, while at larger shaking diameters, i.e. $\tilde {d}_s=0.13$ and $0.20$, a DC wave first emerges on the top of a SC dynamics and eventually a DC wave breaking occurs at larger frequencies, a jump-down transition from DC to SC takes place with increasing $\varOmega$ at lower shaking diameters, i.e. $\tilde {d}_s=0.07$ and $0.10$ for $D=0.144\,\text {m}$. This well-known hysteretic behaviour can be reasonably ascribed to the viscous dissipation of the system. For instance, at sufficiently small shaking diameters, e.g. $\tilde {d}_s\approx 0.02$ (see figure 2), the DC dynamics does not manifest at all, as the energy pumped into the system by the external forcing is likely not sufficient to dominate over the system viscous dissipation, whose effect also depends on the container diameter $D$. Indeed, figure 7 clearly shows that larger diameters, i.e. $D=0.287\,\text {m}$, generate less dissipation. It follows that for larger $D$, on increasing the driving frequency at a fixed shaking diameter, e.g. $\tilde {d}_s=0.10$, the free surface is more likely to undergo wave breaking, rather than a jump-down transition (see figures 6b and 7c). Obviously, the inviscid model employed here is not capable of predicting the so-called jump-down frequency. In Appendix A, a heuristic viscous damping model is introduced to tentatively overcome these limitations.

Finally, we note that, for frequency moderately far from the super-harmonic resonance, the agreement of the WNL model with experiments and with the linear solution, which behaves well away from $\varOmega \approx \omega _{11}$, progressively deteriorates. This is particularly evident on the lower stable branch and can be ascribed to the fact that the asymptotic model is essentially formalized for a fixed driving frequency, i.e. $\varOmega \approx \omega _{21}/2$. The second-order correction (proportional to $\varLambda F$ and discussed in § 4.2.1) to the leading-order particular solution seems to be sufficient to guarantee a fairly good agreement of the upper stable branch in a relatively wide range of frequency $\varOmega <\omega _{21}/2$. However, for the lower stable branch, i.e. $\varOmega >\omega _{21}/2$, the agreement with non-breaking-wave experiments sufficiently far from resonance is still relatively poor, despite the fact that these measurements essentially follow the linear prediction.

4.2.3. The $\epsilon ^3$-order correction to the leading-order SC particular solution

The last consideration in § 4.2.2 suggests that the leading-order SC solution with its second-order correction is only accurate in a limited range of frequency and higher-order terms should be accounted for in order to better retrieve the linear prediction far from resonance. In this regard, as the present asymptotic expansion is pursued up to the $\epsilon ^3$ order, we note that the third-order forcing contains a term, namely

(4.26)\begin{equation} \boldsymbol{\mathcal{F}}_3^{\varLambda\varLambda F}={-}\text{i}\varLambda^2F\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_2^{\varLambda F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)})\exp({\text{i}\varLambda T_1})+{\rm c.c.} , \end{equation}

generated by the derivative of the second-order solution with respect to the slow time scale $T_1$. This term is not resonating in $\exp ({\text {i}(\omega _{2n}-2\theta )})$ and therefore it can be gathered together with the other third-order non-resonating terms (NRT) in (4.19), which in asymptotic models are typically ignored, unless one aims to proceed to higher orders. Nevertheless, such a forcing term produces a system response:

(4.27)\begin{equation} \boldsymbol{q}_3^{\varLambda\varLambda F}=\varLambda^2F\hat{\boldsymbol{q}}_3^{\varLambda \varLambda F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)})\exp({\text{i}\varLambda T_1})+{\rm c.c.} , \end{equation}

which precisely represents the third-order frequency correction of the leading-order particular solution and acts similarly to the second-order frequency correction $\boldsymbol {q}_2^{\varLambda F}=\varLambda F\hat {\boldsymbol {q}}_2^{\varLambda F}\exp ({\text {i}((\omega _{2n}/2)t-\theta )})\exp ({\text {i}\varLambda T_1})+{\rm c.c.}$ discussed in § 4.2.1.

A better intuition about why this third-order correction should improve the prediction farther away from the super-harmonic resonance is given in the following. As anticipated in § 4.2.1, the first-order particular solution simply represents the linear system response to the external forcing $\epsilon F\cos {(\varOmega t-\theta )}$. In general, the resulting amplitude is $\propto \epsilon F/(\varOmega ^2-\omega _{1n}^2)$, but in our asymptotic analysis the leading-order forcing frequency is frozen to $\omega _{2n}/2$, which implies that the first-order particular solution has an amplitude fixed to $\epsilon F/(\omega _{2n}^2 /4-\omega _{1n}^2)$, which does not account for the frequency detuning. If one replaces $\varOmega =\omega _{2n}/2+\epsilon \varLambda$ in $\epsilon F/(\varOmega ^2-\omega _{1n}^2)$ and takes its Taylor expansion, the $\epsilon ^2$-order term is $\propto \epsilon ^2 F\varLambda$, while the $\epsilon ^3$-order term is $\propto \epsilon ^3F\varLambda ^2$. It naturally follows that accounting for this third-order correction (4.27) leads to a more accurate description of the linear response to the external forcing and, therefore, it should give a better model prediction farther away from resonance, where the DC amplitude $|B|\approx 0$ and the non-breaking-wave experiments follow the linear theory.

To conclude, taking the total solution as

(4.28)\begin{equation} \boldsymbol{q}_{DC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2+\boldsymbol{q}_3^{\varLambda \varLambda F} \end{equation}

is expected to leave the amplitude saturation prediction for the DC wave, $|B|$, unaltered, as it does not contribute to the amplitude equation solution, and, simultaneously, to better describe the SC swirling farther away from $\varOmega \approx \omega _{2n}/2$ (at least where no breaking waves occur).

The influence of the aforementioned $\epsilon ^3$-order corrections on the prediction for DC waves is shown as violet lines in figures 6 and 7, where it can be seen that accounting for the additional term allows one to eventually close the gap with the experiments even farther away from resonance, hence ensuring to the WNL model a wider frequency range of validity in all cases examined.

4.2.4. Experiments versus WNL prediction: free-surface reconstruction

In figure 8, the WNL model prediction for the DC waves is compared with the straightforward asymptotic prediction discussed in § 3 and the experimental measurements for DC waves from Reclari (2013) and Reclari et al. (2014). The direct quantitative comparison is here outlined in terms of dimensionless and phase-averaged wave height measured at the sidewall, $\tilde {\delta }(\theta )$.

Figure 8. (a,c,e,g) Comparison of the dimensionless and phase-averaged wave height measured at the wall $\tilde {\delta }(\theta,{\rm \pi} /\varOmega )$ (black circles) with the straightforward asymptotic solution rebuilt via (3.14) (grey solid line) and the WNL solution for the DC wave (4.25). Panels correspond to $\tilde {H}=0.52$, $\tilde {d}_s=0.11$ and $D=0.144\,\text {m}$. The experimental measurements, here shown as black circles, are available in Reclari (2013), except for those of (e), which are provided in Reclari et al. (2014). Note that (c) the nonlinear prediction has a very large amplitude. (b,d, f,h) Corresponding three-dimensional free-surface deformation, $\eta (r,\theta,{\rm \pi} /\varOmega )$, reconstructed via (4.25). The transition from SC to DC swirling via hardening nonlinearity is clearly visible moving from top to bottom, i.e. for increasing frequency.

We observe that, if at $\varOmega /\omega _{21}=0.490$ both models match satisfactorily the experimental points, as soon as $\varOmega /\omega _{21}=0.5$ is approached, the straightforward asymptotic solution diverges due to the resonant (second-order) super-harmonic term, while the WNL solution predicts correctly the finite-amplitude saturation and the emergence of a DC wave on the top of a SC one. The WNL model for DC waves remains in fairly good agreement even at larger driving frequency, although the increasing phase asymmetry between the two local peaks at $\theta ={\rm \pi} /2$ and $3{\rm \pi} /2$ is not retrieved by the present inviscid asymptotic analysis, where secondary effects, e.g. the phase shift induced by viscous dissipation and influence of other higher modes, as well as stronger nonlinear effects for increasing wave amplitudes are overlooked.

For completeness, the three-dimensional free surface, $\eta (r,\theta,{\rm \pi} /\varOmega )$, is reconstructed through (4.25) and shown in figure 8(b,d,f,h), where, for increasing shaking frequencies, the nonlinear transition from nearly SC to DC swirling is evident.

4.2.5. The Helmholtz–Duffing oscillator analogy

While the Duffing equation is known to capture period-3 and period-1/3 dynamics arising from the cubic nonlinearity (Jordan & Smith 1999; Kalmár-Nagy & Balachandran 2011), as those observed by Bäuerlein & Avila (2021) and occasionally by Reclari et al. (2014), it cannot predict the period-halving (the system responds at a frequency which is twice that of the external forcing) dynamics associated with the super-harmonic resonance investigated in this paper. Therefore, in connection with § 4.1.3, here we aim to identify the simplest possible mechanical oscillator that could mimic, at least from a qualitative perspective, the period-1/2 dynamics studied in this work.

The WNL analysis as well as the straightforward asymptotic model highlighted the crucial role of quadratic nonlinearities emerging at second order and from which the DC dynamics stems. At the same time, the WNL model revealed that second-order terms only are not sufficient to capture all the dynamics features owing to the lack of restoring terms and, therefore, cubic nonlinearities must be retained. These considerations suggest that the DC dynamics could be tentatively described by a driven oscillator with both quadratic (asymmetric) and cubic (symmetric) nonlinear terms, i.e.

(4.29)\begin{equation} \ddot{x}+2\sigma\dot{x}+ x+c_2 x^2 + c_3 x^3 = p\cos{\varOmega t}. \end{equation}

Equation (4.29), also commonly known as the Helmholtz–Duffing equation, has wide applications in engineering problems, such as those related to beams, plates and shells subjected to an initial static curvature (Mirzabeigy, Yazdi & Nasehi 2014; Askari et al. 2011), whose governing equations are reduced to a second-order nonlinear ordinary equation with quadratic and cubic nonlinear terms (Ke, Yang & Kitipornchai 2010; Alijani, Bakhtiari-Nejad & Amabili 2011; Fallah & Aghdam 2011).

Among the diverse asymptotic solutions of (4.29) in different limits (Rega 1995; Benedettini & Rega 1989; Kovacic & Brennan 2011), the most relevant to our work is that of Benedettini & Rega (1989). Within the context of planar nonlinear response of suspended elastic cables to an external excitation, they derived an amplitude equation which is concerned with the first or fundamental super-harmonic excitation, i.e. $\varOmega \approx 1/2$, of (4.29). Their WNL approach is detailed in Appendix C, with the additional assumption of vanishing damping $\sigma =0$. Assuming $2\varOmega =1+\lambda =1+\epsilon \hat {\lambda }$, small nonlinearities, $c_2=\epsilon \widehat {c_2}$ and $c_3=\epsilon ^2\hat {c}_3$, and introducing two slow time scales, one obtains

(4.30)\begin{equation} {\rm d}D/{\rm d}t={-}\text{i}\,({2\lambda} +c_5f^2 )D +\text{i}\,\left(1-{\lambda}\right)c_2f^2/2-\text{i}\,4c_4|D|^2D, \end{equation}

with $C=D{\rm e}^{\text {i}{2\lambda } t}$ and with the auxiliary coefficients $c_4$ and $c_5$ (both functions of $c_2$ and $c_3$) defined in Benedettini & Rega (1989). By comparing term by term, the analogy with (4.22) is evident.

To conclude, although the DC dynamics examined in this paper is intrinsically related to the simultaneous interplay of multiple waves, thus making particularly challenging an accurate one-to-one quantitative comparison with a single-degree-of-freedom mechanical model, (4.30) seems to suggest that the actual inviscid sloshing dynamics in the DC regime may be, at least qualitatively, described by the undamped Helmholtz–Duffing equation (4.29) driven super-harmonically.

5. Conclusion

With regards to orbital-shaken cylindrical containers and, specifically, to the careful experimental campaign reported in Reclari (2013) and Reclari et al. (2014), a WNL analysis via a multiple-time-scale method was formalized in § 4 in order to investigate diverse features of the steady-state free-surface dynamics and, particularly, the DC wave dynamics pertaining at half the frequency of the first $m=2$ natural mode.

After having discussed the substantial limitations of the straightforward expansion procedure proposed by Reclari et al. (2014) and summarized in § 3, the WNL analysis was first formulated under the most common condition of pure harmonic resonance. Despite the inviscid assumption, the WNL analysis developed for the SC wave dynamics was shown to be in fairly good agreement with all the experimental measurements. In fact, the present model correctly describes the close-to-resonance hardening nonlinear behaviour experimentally observed. The agreement remains sufficiently accurate until the free surface eventually breaks and a transition to a fully nonlinear regime occurs.

It is well assessed in the literature that the close-to-harmonic-resonance sloshing dynamics can be modelled (from both qualitative and quantitative perspectives; Bäuerlein & Avila 2021) by a single-degree-of-freedom system with a cubic nonlinearity and driven harmonically, i.e. by the famous Duffing oscillator, as rigorously proved for a two-dimensional rectangular container laterally excited (Ockendon & Ockendon 1973). Without surprise, this was shown to hold for the case of orbital-shaken cylindrical containers as well.

The WNL analysis was then extended to the more complex case of a DC swirling. The overall agreement with experiments and, especially, the improvements with respect to the simple straightforward asymptotic model are remarkable in all cases considered, although the slight asymmetry observed in the reconstruction of the periodic free-surface dynamics at the sidewall was not retrieved in the present model.

To the knowledge of the authors, a formal amplitude equation describing the super-harmonic DC sloshing dynamics in orbital-shaken containers and coupled with a thorough experimental validation has not been reported in the literature yet, hence representing the most significant finding of this work.

Lastly, by analogy with the close-to-harmonic-resonance dynamics for SC waves, for which the Duffing oscillator represents a suitable mechanical analogy, a one-degree-of-freedom mechanical oscillator having both quadratic and cubic nonlinear terms, commonly referred to as a Helmholtz–Duffing oscillator, driven super-harmonically, was tentatively identified as the simplest possible mechanical system that could mimic, at least qualitatively, the super-harmonic DC sloshing dynamics investigated in this paper. The Helmholtz–Duffing equation was largely adopted in the last few decades within the context of structural analysis, i.e. beams, plates and shells subjected to an initial static curvature as well as suspended elastic cables (Nayfeh 1984; Benedettini & Rega 1989), and it was here proposed as a direct mechanical analogy with the present orbital sloshing system.

The main limitation of the models derived in this work is intrinsic to the fundamental assumption of an inviscid fluid. This precludes correctly accounting for the jump-down transition experimentally observed for DC waves at low shaking amplitudes and, therefore, for an accurate estimation of the maximum amplitude response when such a transition occurs. Furthermore, in the absence of viscous boundary layers, the WNL time- and azimuthal-averaged mean flow reduces to a free-surface deformation only. This is in stark contrast with the existence of the so-called Eulerian mean flow (van den Bremer & Breivik 2018), also known as viscous streaming flow, typically observed in experiments (Bouvard et al. 2017). Therefore, the present work overlooks one of the essential points of interest in applications of orbital shaking. The mean flow, which contributes to an efficient mixing, is not captured.

The extension of the asymptotic models developed in this work to a viscous analysis is desirable, as it would enable one to predict quantitatively these secondary but fundamental effects for both cases of harmonic and super-harmonic resonances. However, it presently hinges on the subtle modelling of the moving contact line condition.