3.1. Assessing the quasi-two-dimensionality of the flow

The geometry of our tank suppresses swirling motions and to the naked eye the motion of the sloshing water surface appears quasi-two-dimensional for various driving parameters and sloshing states. In order to quantify this, we performed stereoscopic PIV measurements at the highest investigated sloshing amplitude (in which the liquid rises steeply at the tank wall). The snapshots of the velocity field shown in figure 4(a) were recorded when the surface elevation was minimal (and the flow velocities were maximal). The in-plane velocities ($v_x, v_y$) reflect the sloshing motion, while the velocities in $z$-direction are negligible (and remain so throughout the cycle, see the time series of figure 4b). Note that only every eighth vector obtained from the PIV measurements is plotted here for clarity. In the time series of figure 4(b) the deviations of the $v_z$-component reflect the vibrations of the cameras itself (see § 2.3) and are very small, confirming that the flow is quasi-two-dimensional. For the following analysis of time series (in figure 10 and 11) we applied a low-pass filter to remove these vibrations from the signal. In figure 4(c) it is shown that the in-plane velocities decrease exponentially with the distance from the surface, as it is expected for deep water surface waves. The exponential decrease was observed at all positions.

3.2. Characterising the sloshing amplitude with the liquid's centre of mass

The most common method used in the literature to characterise the sloshing amplitude is to measure the maximum wave elevation, typically at a fixed position in the container (Dodge Reference Dodge2000; Ibrahim Reference Ibrahim2005; Faltinsen Reference Faltinsen2017). This is a direct and simple method and the data can be retrieved from wave sensors or in situ with a ruler. Using our image analysis method, we measure the surface displacement $\zeta$ with respect to the resting liquid height at the sidewall. The amplitudes of the wave crests (maxima) and troughs (minima) are shown in figure 5(a) as a function of the excitation frequency for $A=0.64\,\%$. Far from the resonance (maximum response), the amplitudes of the waves and crests are similar, but close to the resonance the crests are very steep and the troughs are broad (Dodge Reference Dodge2000; Ibrahim Reference Ibrahim2005), leading to distinctly different amplitudes (as it is illustrated by the arrows in the snapshot in (c)). In many studies, it is not explicitly stated whether crests or troughs are used to characterise the sloshing amplitude, which hinders quantitative comparisons. Clearly, the surface elevation at a fixed location does not provide global information on the sloshing state, and as it will be demonstrated in the next section, this can cause problems when measuring the viscous damping coefficient.

Figure 5. Response curves when starting the experiment from rest at $A=0.64\,\%$. (a) The surface displacements of wave maxima and minima (crests and troughs) lead to different response curves of the sloshing amplitude $\zeta$. (b) The amplitude $\hat X$ of the horizontal oscillation of the liquid's centre of mass provides a global measure of the response curve and allows a precise determination of the sloshing phase. The solid line (grey) corresponds to a harmonic oscillator (without fitting parameters) and the red line to a Duffing oscillator (one fit parameter, unstable solutions are denoted with dotted segments). (c) Flow snapshot (with the free surface shown as a red line) illustrating the typical asymmetry in the surface displacement close to resonance (at $\varOmega = 0.953$). This asymmetry causes the different response curves shown in (a). The liquid's centre of mass is displayed as a red dot.

In engineering applications, the motion of the centre of the liquid's mass is often the main quantity of interest (e.g. for the propellant tank of a satellite). The identification of the water surface combined with the quasi-two-dimensionality of the flows allows us to determine directly the centre of the liquid's mass (see the red dot in figure 5c) as a function of time, $x(t)$, which we use hereafter as global measure of the sloshing amplitude. More specifically, we use the amplitude of the horizontal periodic displacement of the liquid's centre of mass, termed $\hat x$ (or $\hat X = \hat x/w$ for the dimensionless amplitude) hereafter. The corresponding response curve is shown in (b). Being based on the analysis of the whole surface shape, this method is robust against small image evaluation errors and is well suited to quantitatively compare different flow states.

3.3. Determining the damping rate

The viscous damping rate $\delta _1$ is an important parameter required in models of sloshing, including the Duffing oscillator and multimodal models. In this section, it is determined experimentally and analytically for our system. In experiments, the viscous damping rate $\delta _1$ can be extracted from the system's natural response. After turning off the driving, the sloshing amplitude decays over time to zero with decay rate $\delta _1$. This is determined by fitting an exponential function to the envelope of the time series of the sloshing amplitude. In aerospace engineering it is common to use the ‘free surface displacement at some convenient location’ (Dodge Reference Dodge2000) to obtain the damping coefficient. In figure 6(a) it is shown that this method leads to different damping coefficients for the crests and the troughs. By contrast, the motion of the liquid's centre of mass allows an unambiguous determination of the viscous damping coefficient for each time series (see figure 6b). Repeating such measurement for a range of excitation frequencies revealed only small deviations, $\delta _1 = [0.065\pm 0.015]\ \mathrm {s}^{-1}$ (corresponding to dimensionless damping ratio $\gamma _1=\delta _1/\omega _1=8.4 \times 10^{-3}$). For very small frequencies (indicated by $\circ$ (red) in figure 7a) the sloshing amplitude was too small (surface wave height $< 6\ \mathrm {mm}$) to accurately fit a decaying function over the signal. As a result, these values were not considered to compute the averaged damping rate $\delta _1$. Keulegan (Reference Keulegan1959) derived an analytical estimate of the damping rate for liquid sloshing in rectangular containers. In his model, the damping rate is obtained from the balance of potential energy, kinetic energy and energy dissipation caused by a laminar viscous boundary layer on the wetted tank walls and bottom. Keulegan's prediction reads (Faltinsen & Timokha Reference Faltinsen and Timokha2009),

(3.1)\begin{equation} \gamma_n =\frac{\delta_n} {\omega_n}= \sqrt{\frac{\nu}{2 w^{2} \omega_n}} \left[\left(\frac{w}{l}\right) \left( 1 + 2{\rm \pi} \frac{\displaystyle n \left(\frac 12 - \frac{h}{l}\right)}{\displaystyle \sinh\left[2{\rm \pi} n \frac{h}{l} \right]} \right) +1 \right], \end{equation}

where $\gamma _n$ is the dimensionless damping ratio of the natural mode $n$. In our decay experiments the first natural mode $n=1$ dominates. Plugging the kinematic viscosity $\nu = 1 \times 10^{-6}\ \mathrm {m^{2} /s}$ for water at room temperature and our experimental geometry ($h=0.4\ \mathrm {m}$, $w=0.5\ \mathrm {m}$, $l=0.05\ \mathrm {m}$) into (3.1) gives $\gamma _1 = 5.57 \times 10^{-3}$ ($\delta _1=0.0434 \ \mathrm {s}^{-1}$), which is approximately $1.5$ times smaller than the experimentally determined value, $\gamma _1=8.4 \times 10^{-3}$. This difference is unsurprising, as (3.1) only models the effect of a linear boundary layer and does not account for damping in the fluid bulk, turbulent energy dissipation, surface tension or free-surface effects like wave breaking (Faltinsen & Timokha Reference Faltinsen and Timokha2009). Moreover, nonlinear effects may be relevant in the early stages of the decay, making it eventually necessary to consider also the damping of higher natural modes. As a consequence, the theory only gives a lower bound for the actual damping. Similar differences were observed in past studies (Ikeda et al. Reference Ikeda, Ibrahim, Harata and Kuriyama2012; Faltinsen & Timokha Reference Faltinsen and Timokha2017). In our case it is noteworthy that Keulegan's prediction agrees with our measurements at low sloshing amplitude (corresponding to a low excitation frequency in figure 7a), where the theory is expected to be most accurate. However, we stress that damping rates at low sloshing amplitudes could not be measured accurately. In the following, we use the experimentally determined damping ratio $\gamma _1=8.4 \times 10^{-3}$ for comparison with sloshing models.

Figure 6. Dependence of the viscous damping rate on the observable. At $t=0$ the driving was turned off to measure the decay rate. The steady state ($t\le 0$) is a wave-breaking motion, which fades into a planar surface wave during the decay ($t>0$). (a) Time series of the dimensionless surface displacement, $\zeta /w$, recorded at the left tank wall (where the sloshing amplitude is maximal) for $A = 0.64\,\%$, $\varOmega = 0.999$. The asymmetry of the oscillation amplitude reflects the asymmetric wave shape (with a higher wave crest and a shallow wide trough), leading to different positive and negative displacements ($C_a = 0.153w$, $C_b = 0.126w$) and hence different decay rates ($\delta _a = 0.069\ \mathrm {s}^{-1}$, $\delta _b = 0.060\ \mathrm {s}^{-1}$). (b) The same time series as in (a) but for the lateral motion of the centre of mass. The displacement of the centre of mass is left–right symmetric and an unambiguous damping coefficient is obtained ($C = 0.035w$, $\delta _1 = 0.064\ \mathrm {s}^{-1}$).

Figure 7. Measured viscous damping rate and natural frequency. (a) Experimentally determined values of the viscous damping rate $\delta _1$ for different excitation frequencies and $A=0.64\,\%$. Their average $\delta _1= 0.065\,1/s$ is shown as a red line. (b) Typical frequency spectrum of a decaying oscillation of the liquid's centre of mass with a peak at $\omega _1$, showing that the decaying motion oscillates with the natural frequency. The spectrum stems from the time series of the decay shown in figure 6(b).

3.4. Determining the natural frequency

The natural frequency of the sloshing can be calculated analytically with potential theory (2.2) to $\omega _1= 7.800\ \mathrm {s}^{-1}$. Experimentally, we verified its value with the same measurements that were also used to obtain the viscous damping rate $\delta _1$. After stopping the forcing, the system oscillates with its natural frequency, which is seen in the frequency spectrum in figure 7(b). The averaged value of the measurements ($\omega _{1, exp}=[7.86\pm 0.10]\ \mathrm {s}^{-1}$) agrees well with the prediction. This confirms that the measured damping rate can be attributed to the first natural mode.

4.1. Fitting the nonlinearity constant

The nonlinearity constant $\epsilon$ is the only parameter in (4.7) which can neither be measured, nor estimated a priori with potential theory. The red line in figure 5(b) shows a least-square-residuals fit of (4.7) to the horizontal displacement of the liquid's centre of mass, performed simultaneously for all data sets measured in this work (see § 7). This demonstrates that, even at the large amplitude investigated here, the Duffing equation captures the measured response very well (with a single fit parameter, $\epsilon = -1.44 \times 10^{4} \ \mathrm {m}^{-2}\ \textrm {s}^{-2}$). Note that the harmonic oscillator model ($\epsilon =0$), which is entirely free of fitting parameters, is able to very precisely capture the system's response far away from the resonance, see the grey line in figure 5(b).

In order to ease future comparisons with numerical simulations and experiments with other dimensions, all parameters are shown hereafter in dimensionless form. With the dimensionless horizontal displacement of the liquid's centre of mass, $X=x/w$, the dimensionless Duffing equation reads

(4.9)\begin{equation} \ddot X + 2\gamma_1 \dot X + X + \varepsilon X^{3} = c_1 A \varOmega^{2} \cos \varOmega t, \end{equation}

where time has been rescaled with the inverse of the natural frequency $\omega _1$, $\gamma_{1} = \delta _1 / \omega _1 = 8.4 \times 10^{-3}$ is the dimensionless damping coefficient and $\varepsilon =(w/\omega _1)^{2}\epsilon =-59.2$. A comparison between the dimensional and dimensionless parameters is shown in table 1.

Table 1. Experimentally determined parameters of the Duffing equation in dimensional (4.5) and dimensionless (4.9) forms. For the measurement procedure, see § 4.

5.1. Model equations for sloshing in a rectangular tank

Faltinsen & Timokha (Reference Faltinsen and Timokha2009) give a solution of the potential-flow equations for two-dimensional inviscid sloshing in a rectangular container using a Fourier expansion. The elevation of the free surface as a function of the horizontal container coordinate, $\xi \in [-w/2, w/2]$, and of time, $t$, can be expressed as

(5.1)\begin{equation} \zeta(\xi,t) = w \sum_{n=1}^{\infty} \beta_n(t) f_n(\xi), \end{equation}

where

(5.2)\begin{equation} f_n(\xi) = \cos({\rm \pi} n (\xi+ \tfrac 12 w)/w ), \end{equation}

is the dimensionless waveform of the surface for mode $n$ and $\beta _n$ describes its temporal evolution. When this expression is inserted in the potential-flow equations, these reduce to an infinite system of (coupled) nonlinear differential equations for the time-dependent coefficients $\beta _n$. By truncating this expansion, multimodal models of sloshing of desired order can be obtained. Faltinsen & Timokha (Reference Faltinsen and Timokha2009) provide equations for $\beta _n$ for the first three modes ($n=1,2,3$), which constitute the lowest-order, consistent nonlinear multimodal model of sloshing in a rectangular container (Moiseev Reference Moiseev1958). The analytic solution to these equations is given in chapter 8.3 of Faltinsen & Timokha (Reference Faltinsen and Timokha2009) and reads

(5.3)\begin{equation} \left. \begin{array}{c@{}} \beta_1(t) = \hat Y \cos(\omega t) + \tilde{n}_1 \hat Y^{3} \cos(3\omega t) + O(\hat Y^{5}), \\ \beta_2(t) = \hat Y^{2} (w_0 + h_0 \cos(2\omega t)) + O(\hat Y^{4}), \\ \beta_3(t) = \cos(\omega t) [\tilde{N}_1 \hat Y^{3} - A P_3 / (1-(\omega_3 / \omega)^{2}) ] + \tilde{N}_2 \hat Y^{3} \cos(3\omega t) + O(\hat Y^{5}), \end{array}\right\} \end{equation}

where $\hat Y$ is a dimensionless amplitude, $\omega _n$ are the natural frequencies of the three modes (see (2.2)). Hence, the time-dependent coefficients $\beta _n$ of the solution (up to third order in amplitude $\hat Y$) consist of harmonics of the driving frequency $\omega$, whose magnitude depends on the dimensionless amplitude $\hat Y$, and the following parameters:

(5.4)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \tilde{n}_1 ={-} \dfrac{d_2 + h_0(3d_1 + 4d_3)}{2(9- (\omega_1 / \omega)^{2})}, \\ \displaystyle \tilde{N}_1 = \dfrac{-3q_2 + q_4 + 2h_0({-}q_1 - 4q_3 + 2q_5) - 4q_1 w_0}{4(1 - (\omega_3 / \omega)^{2})}, \\ \displaystyle \tilde{N}_2 ={-} \dfrac{q_2 + q_4 + 2h_0(q_1 + 4q_3 + 2q_5)}{4 (9 - (\omega_3 / \omega)^{2})}, \end{array}\right\} \end{equation}

and

(5.5a,b)\begin{equation} w_0 = \frac{d_4 - d_5}{2 (\omega_2 / \omega)^{2}}, \quad h_0 = \frac{d_5 + d_4}{2((\omega_2 / \omega)^{2} - 4)}. \end{equation}

The coefficients $d_i$ and $q_i$ depend on the ratio of filling level to container width, $h/w$, and are defined as,

(5.6)\begin{equation} \left.\begin{array}{c@{}} \displaystyle d_1 = 2 \dfrac{E_0}{E_1} + E_1, \quad d_2 = 2 E_0 \left({-}1 + \dfrac{4 E_0}{E_1 E_2} \right), \\ \displaystyle d_3 ={-}2 \dfrac{E_0}{E_2} + E_1, \quad d_4 ={-}4 \dfrac{E_0}{E_1} + 2 E_2, \\ \displaystyle d_5 = E_2 - 2 \dfrac{E_0 E_2}{E_1^{2}} - 4 \dfrac{E_0}{E_1}, \end{array}\right\} \end{equation}

and

(5.7)\begin{equation} \left.\begin{array}{c@{}} \displaystyle q_1 = 3 E_3 - 6 \dfrac{E_0}{E_1}, \quad q_2 ={-}3 E_0 - 9 \dfrac{E_0 E_3}{E_1} + 24 \dfrac{E_0^{2}}{E_1 E_2}, \\ \displaystyle q_3 = 3 E_3 - 6 \dfrac{E_0}{E_2}, \quad q_4 ={-}6 E_0 - 24 \dfrac{E_0 E_3}{E_1} + 48 \dfrac{E_0^{2}}{E_1 E_2} + 24 \dfrac{E_0^{2} E_3}{E_1^{2} E_2}, \\ \displaystyle q_5 = 6\left(\dfrac 12 E_3 - \dfrac{E_0}{E_1} - \dfrac{E_0 E_3}{E_1 E_2} - \dfrac{E_0}{E_2}\right), \end{array}\right\} \end{equation}

respectively, with

(5.8)\begin{equation} \left.\begin{array}{c@{}} \displaystyle E_0 = \dfrac{1}{8} {\rm \pi}^{2}, \quad E_n = \dfrac 12 {\rm \pi}\tanh \left[{\rm \pi} n \dfrac{h}{w} \right], \quad n \geq 1 \\ \displaystyle P_n = \dfrac{2}{n {\rm \pi}} \tanh \left[{\rm \pi} n \dfrac{h}{w}\right] (({-}1)^{n} - 1). \end{array}\right\} \end{equation}

Under the assumption of periodic (steady-state) oscillations and by further considering linear damping added to the equation for the first mode, Faltinsen & Timokha (Reference Faltinsen and Timokha2017) obtained an analytic expression for the dimensionless amplitude parameter

(5.9)\begin{equation} \hat Y = \frac{P_1 A}{\sqrt{[ (\omega_1/\omega)^{2} - 1 + M_1 \hat Y^{2} ] + \gamma^{2}}}, \end{equation}

where

(5.10)\begin{equation} M_1 = d_1 ({-}w_0 + \tfrac 12 h_0) - \tfrac 12 d_2 - 2 d_3 h_0. \end{equation}

The phase lag is

(5.11)\begin{equation} \cos\phi = \frac{\hat Y [(\omega_1/\omega)^{2} - 1 + M_1\,\hat Y^{2} ]}{P_1 A}. \end{equation}

Equations (5.9) and (5.11) are used to generate the response curves of the multimodal model for our system parameters.

5.2. Similarities and differences between the multimodal model and the Duffing equation

Equations (5.9) and (5.11) for the amplitude $\hat Y$ are very similar in structure to the solution of the Duffing equation for the amplitude of the liquid's centre of mass $\hat X$ (shown in this paper only in dimensional form, i.e. (4.7) and (4.8) for $\hat x$). The square of the amplitude parameter $\hat Y$ is multiplied by a nonlinearity parameter $M_1$ akin to the Duffing parameter $\epsilon$. The response curves of the multimodal model are thus qualitatively similar to those of the Duffing oscillator, as will be shown later in § 7. For example, as the driving amplitude increases, the resonance peak bends and a hysteretic region with two stable solutions and an unstable solution develops. Depending on the sign of $M_1$ this resonance can be softening (bending towards smaller frequencies) or hardening (towards higher frequencies). More specifically, at a critical filling of $h_\text {crit}/w = 0.3368$ the coefficient becomes zero (Ockendon & Ockendon Reference Ockendon and Ockendon1973), with softening (hardening) predicted for larger (smaller) $h/w$. In our case $h/w=0.8$ and accordingly our system exhibits a softening response. A clear advantage of the multimodal model over the Duffing equation is that $M_1$ can be computed a priori, whereas $\epsilon$ must be obtained with a fit to our experimental data. Note also that $M_1$ is a monotonically decreasing function of the driving frequency $\omega$. In the next section we compare experimentally measured surface shapes to the predictions of the mulitmodal model of Faltinsen & Timokha (Reference Faltinsen and Timokha2017).

5.3. Comparison of the predicted and measured surface wave shapes

The multimodal model can be used to predict the temporal evolution of the surface shape during an oscillation period (under steady-state conditions). For our geometry and selected driving parameters, the amplitude parameter $\hat Y$ is computed according to (5.9), the time-dependent functions $\beta _n$ are evaluated according to (5.3) and then inserted into (5.1), yielding a modal prediction of the wave elevation $\zeta$. Three examples of calculated waveforms are shown in the upper row in figure 8 for varying excitation frequencies. For $\varOmega = 0.950$, shown in (a), the amplitude parameter is $\hat Y=0.084$ and the wave shape is clearly dominated by the first mode (shown as dotted line). As the second mode is added, the prediction (dash-dotted line) becomes slightly better, but it does not improve much when the third mode is also considered (solid line). For $\varOmega = 0.962$ in (c), the amplitude parameter is much larger, $\hat Y=0.233$, and as a consequence the contributions of the second and third modes become more significant, which results in a steeper wave shape, a higher wave crest and a shift of the nodal point away from the centre. These features are in qualitative agreement with our experimental measurements, however, three distinct differences are identified. In the experiments the rising of the liquid at the wall is substantially steeper, the trough is slightly more shallow and a local maximum of the surface inside the trough is not observed. For larger sloshing amplitude, $\varOmega = 0.909$ and $\hat Y = 0.377$, these differences are strongly exacerbated (as shown in $(\textit {e})$), leading to very different crest maxima and qualitatively different wave shapes (with less maxima/minima in the experiments). This suggests that additional (higher-order) modes or stronger damping (especially for the modes $n=2$ and $n=3$, where damping is neglected) may be necessary to accurately capture the behaviour of the system.

Figure 8. Comparison of the free-surface elevation from the multimodal model in (5.1) to the experimentally obtained one for varying driving frequency $\varOmega$ and constant amplitude $A=0.64\,\%$. The black dashed lines illustrate the experimentally obtained surface and the coloured lines the ones of the multimodal model for increasing number of modes as indicated in the legend in a. The red solid line denotes e.g. the predicted surface elevation based on three modes ($N=3$). Snapshots of the surface elevation over the container width $\xi$ are shown in the upper panel and the corresponding time series at the left wall in the panel below. For a better visual comparison, the phases of the oscillations between experiment and model are optimally aligned. The driving parameters are in (a,b) $\varOmega = 0.950$, $\hat Y = 0.084$, in (c,d) $\varOmega = 0.962$, $\hat Y = 0.233$ and in (e,f) $\varOmega = 0.909$ and $\hat Y = 0.377$.

In the literature, the complete surface shape has rarely been quantitatively analysed experimentally, which makes it difficult to improve and validate models. Typically, the time series of the elevation at a fixed location, e.g. at the tank wall, is used for comparison. For completeness, the corresponding time series for the three cases above are shown in figure 8(b,d,f). For $\varOmega = 0.950$ in (b) and $\varOmega = 0.962$ in (d) an excellent agreement is obtained in both cases despite the marked differences when the whole surface is considered. Hence, our results indicate that single-point measurements might not be sufficient to reliably evaluate model predictions. We conclude that, at low sloshing amplitude, the wave shape of the third-order multimodal model is in excellent agreement with the experiments, whereas at larger amplitudes the wave shape differs substantially, even if the time series at a particular points still yields a good agreement. At larger sloshing amplitudes (closer to the nonlinear resonance and displayed in (f)), the time series and the wave shape deviate more substantially.

6.1. Hysteresis

As customary in the literature (Fultz Reference Fultz1962; Arndt & Dreyer Reference Arndt and Dreyer2008; Konopka et al. Reference Konopka, De Rose, Strauch, Jetzschmann, Darkow and Gerstmann2019), we started our experiments from rest and waited until transients lapsed to obtain a response curve. Subsequently, we determined the amplitude of the sloshing waves by measuring the amplitude of the horizontal periodic displacement of the liquid's centre of mass $\hat X$. As shown in figure 5(b) for $A=0.64\,\%$, this procedure yields a resonance curve with a slight asymmetry (softening), which can also be observed in the measurements of the wave height at the sidewall shown in (a). However, starting the system form rest is insufficient to fully characterise nonlinear resonances (see e.g. Strogatz Reference Strogatz2018). In fact, a different picture emerged when the excitation frequency was changed quasi-statically in a sweep-up, sweep-down procedure (see figure 9a). This allowed it to track stable flow states through the parameter space and revealed the occurrence of a pronounced hysteresis. Compared with the measurements started from rest, the maximum sloshing amplitude occurs at smaller driving frequencies and is substantially higher.

Figure 9. Response curve when performing frequency sweeps at $A=0.64\,\%$. (a) Sloshing amplitudes obtained with quasi-static changes of the excitation frequency during a single run are marked with $\circ$ (blue) for the frequency sweep-up and with $\bullet$ (red) for the frequency sweep-down. The observed response curve ($\times$), when experiments were started from rest (same measurements as in figure 5b) are included for comparison. Arrows indicate the (jump-up/down) transitions between the lower (low sloshing amplitude) and upper (high sloshing amplitude) branches. Hysteresis occurs in between the arrows. The line - $\cdot \ \cdot$ - represents the fitted amplitude curve from the Duffing oscillator, with unstable solutions as dotted segments. (b) Snapshots of different sloshing states recorded at the marked points in (a); i: planar motion, ii: wave breaking, iii: period-three motion. The corresponding movies of the three states (supplementary movies 1, 2, 4 are available at https://doi.org/10.1017/jfm.2021.576) are provided online.

6.2. Classification of flow states

At low excitation frequencies ($\varOmega \lesssim 0.9$) the sloshing state consists of a quasi-planar wave, as exemplified in the snapshots of figure 9(b.i) and figure 10(a). Its linear nature can be observed in the velocity time series in figure 10(e), and more precisely in the frequency spectrum in (i) dominated by the driving frequency and its higher harmonics. At low excitation frequencies, the quasi-planar wave (termed also the lower branch state hereafter) is the only stable flow state and is obtained independently of the initial conditions, i.e. start from rest or sweep-up, sweep-down. Similarly, at large frequencies ($\varOmega \gtrsim 1$) there is only one stable upper branch state, distinguished by its plateau-like wave crest near the wall, as displayed in figure 10(b). The frequency spectrum is similar, except that the second harmonic can hardly be discerned, while the third harmonic is rather strong. For modelling approaches, these frequency spectra provide relevant information. In the multimodal model, the second harmonic, with amplitude $O(\hat Y^{2})$, is included in $\beta _2$ (see (5.3)), whereas this does not appear in the Duffing equation, which features terms of amplitude $O(\hat X)$ and $O(\hat X^{3})$ only. The weakening of the second harmonic for increasing sloshing amplitude, combined with the strengthening of the third harmonic (which appears naturally in the Duffing equation), indicates that a reduction of the dynamics to a spectral submanifold (e.g. represented by the Duffing equation) might be possible (Breunung & Haller Reference Breunung and Haller2018).

Figure 10. Characterisation of the sloshing states at $A=0.64\,\%$. (a)–(d) Exemplary snapshots of the states. The surface is indicated with ——– (red) and the surface height at rest with ······ (orange). (e)–(h) Time series of the horizontal velocity $v_x$ in the bulk of the flow obtained with PIV (recorded at least at $200\ \mathrm {Hz}$). The velocity is spatially averaged over ${\rm \Delta} x=28\ \mathrm {mm}$ and $\varDelta y=28\ \mathrm {mm}$ with the central position being at $x=0$ and $y=0.125h$. (i)–(l) Spectra of the time series in (e)–(h) obtained with fast Fourier transforms. The investigated flow states are: (a,e,i) planar motion $\varOmega =0.944$, (b,f,j) plateau peaks $\varOmega =1$, (c,g,k) wave breaking $\varOmega =0.926$, (d,h,l) peaks $\varOmega =0.934$.

In the following, we start from this plateau-like state (at $\varOmega \sim 1$) and describe the changes in the dynamics as the excitation frequency is decreased on the upper branch. Progressively, the surface plateau becomes broader and extends towards the centre and its leading front becomes increasingly steeper. As the frequency is further decreased, a near-vertical liquid column forms and finally it overturns, causing wave breaking, see the snapshots in figure 9(b.ii) and 10(c). We define the onset of this wave-breaking state, when splashing is first observed. Its frequency spectrum is broader, reflecting the enhanced fluctuations and being suggestive of chaos. For a further reduction of the driving frequency, the wave becomes higher (without breaking) and the liquid rises steeply at the tank wall instead of overturning. We refer to this flow state shown in figure 10(d) as the peak state. Its frequency spectrum is again less broad and qualitatively similar to the plateau state at substantially lower sloshing amplitudes.

For a further decrease of the driving frequency, a remarkable phenomenon is observed: a period-three motion where every third wave peak rises higher at the tank wall. This effect is clearly visible by naked eye (see the sequence of snapshots in figure 11a) and leads to a sub-harmonic frequency $\omega /3$ and additional linear combinations in the frequency spectrum in figure 11(c). The resulting period-three motion is also visualised in the phase portrait in figure 11(d). Interestingly, $3\omega$ and $\omega /3$ are classical secondary frequencies of the Duffing oscillator, stemming from its cubic nonlinearity (Kalmár-Nagy & Balachandran Reference Kalmár-Nagy and Balachandran2011). To our knowledge, these have not been observed experimentally so far, neither in sloshing fluids nor in simple systems mimicking Duffing oscillators. After a sufficient reduction of the driving frequency the period-three motion transitions (at the so-called jump-down transition) to the lower branch with the quasi-planar waves of low amplitude.

Figure 11. Period-three motion at $A=0.64\,\%$ and $\varOmega =0.917$. (a) Snapshots of three consecutive peaks of the oscillation. (b) Time series of the horizontal velocity $v_x$ in the bulk of the flow. The velocity is spatially averaged over ${\rm \Delta} x=28\ \mathrm {mm}$ and ${\rm \Delta} y=28\ \mathrm {mm}$ with the central position being at $x=0$ and $y=0.125h$ obtained with PIV (recorded at $300\ \mathrm {Hz}$). (c) Spectra of the time series in (b) obtained with a fast Fourier transform. (d) Phase portrait of the dynamics. The measurements correspond to 30 periods and are plotted as symbols, but appear as a line due to the low scatter in the data. A movie of the period-three motion (supplementary movie 4) is provided online.

6.3. Co-existence of upper branch states in parameter space

Several repetitions of the above described measurements reveals that the flow transitions are very sensitive to experimental conditions, and that the peak state and the wave-breaking state coexist in a large span of driving frequencies, as shown in the response diagram in figure 12(a). These flow states seem to be marginally stable, we even occasionally observed a change from the peak state to the wave-breaking state without any change in the driving parameters. On the other hand, the plateau state at large driving frequencies ($\varOmega \gtrsim 1$) and the period-three motion at resonance are robust and were always observed. We speculate that further stable and also unstable states (solutions of the Navier–Stokes equations) may exist in phase space. A full understanding of the high-dimensional dynamics observed here poses a challenge to be tackled with the direct numerical computation of periodic solutions of the governing equations (Kawahara, Uhlmann & Van Veen Reference Kawahara, Uhlmann and Van Veen2012).

Figure 12. Nonlinear competition of sloshing states at $A=0.64\,\%$. (a) Measured sloshing amplitudes as function of the driving frequency and (b) corresponding phase lag between sloshing motion and driving. The different symbols represent the observed flow state for each measurement point. The measurements in this figure were obtained with different methods (starting from rest and frequency sweeps with variable frequency increments). Differentiation between states is best seen in the phase difference. In this visualisation lower and upper branches appear exchanged, because of the strongly negative phase lag characteristic of the large-amplitude sloshing branches.

The amplitude of the analytic solutions of the Duffing oscillator (4.7) are plotted as lines in figure 12(a). The measurements in (a) follow closely the stable solutions (solid lines) with the important exception that the jump-down transition (and therefore the response maximum) occurs at substantially higher driving frequencies in the experiment. The jump-down transition obtained from the Duffing oscillator fit is at $\varOmega \approx 0.78$ with an amplitude of $\hat X=9.44\,\%$. In the sloshing experiments such high amplitudes are not reached, suggesting that enhanced dissipation prevents the development of such flow states.

In undamped (inviscid) systems the phase lag between excitation and response is zero. The phase lag of the Duffing oscillator (4.8) is plotted as lines in figure 12(b), together with the experimentally measured values. Note that, only the sloshing amplitude was used to obtain the fitting parameters and not the phase difference. It can be seen that the phase lag of the low-amplitude sloshing with quasi-planar waves (marked as green triangles), as well as of the plateau state, are again in excellent agreement with the Duffing oscillator fit. However, the phase lag of the larger-amplitude sloshing states (peak state, wave breaking, period three) substantially deviates from the Duffing curve. This suggests that the damping is nonlinear for these complex flow states.

The transition from the period-three motion – dominating close to the response maxima – to the quasi-planar waves occurs at approximately $90^{\circ }$. Indeed, the phase lag measurements shown in figure 12(b) evidence that there are two distinct period-three states – both exhibiting the jump-down transition at about $90^{\circ }$.

We conclude that the sloshing amplitude (defined by the liquid's centre of mass) is neither suited to distinguish between flow states, nor does it act as an early indicator for the response maxima – as the jump-down transition appears without warning. The phase lag, on the other hand, allows a better distinction of the flow states and of the deviations from Duffing dynamics as the jump-down transition is approached. It seems that the jump-down transition occurs at a phase lag of $90^{\circ }$, which we further test at different driving amplitudes in the next section.

7.1. Comparison of the experiments to the response curves of a Duffing oscillator

In this section, we quantitatively compare our experimentally obtained response curves for three additional driving amplitudes, $A=0.32\,\%$, $0.17\,\%$ and $0.09\,\%$ to the Duffing oscillator. The only parameter of the Duffing equation in (4.9) that could neither be obtained from potential theory nor measured directly is the nonlinearity constant $\varepsilon =-59.2$, which was freely fitted. The response curves are shown in figure 13, where experimental measurements are represented by symbols and the solid (dotted) lines denote stable (unstable) solutions of the Duffing equation.

Figure 13. Quantitative comparison between sloshing experiments (symbols) and Duffing oscillator (lines) with a single free fitting parameter, $\varepsilon =-59.2$ in (4.9). (a) Sloshing amplitude $\hat X$ exhibiting the characteristic bending of the resonance curve of a softening spring and the development of hysteresis for increasing excitation amplitude $A$. Solid (dotted) lines denote stable (unstable) solutions of the Duffing equation. The measurements (symbols) were obtained by quasi-statically changing the excitation frequency $\varOmega$, whilst keeping the excitation amplitude $A$ fixed. The experimentally determined maximum sloshing amplitudes are indicated by crosses ($\times$). At $A=0.64\,\%$ the upper branch is composed of at least three branches of competing states (see figure 12). (b) Phase lag between sloshing and excitation. The transition from high- to low-amplitude sloshing occurs at a phase difference of $\phi \approx -90^{\circ }$.

At the lowest driving amplitude investigated, quasi-planar waves are found across the whole frequency range. The system behaves here nearly like a harmonic oscillator, but with the maximum sloshing amplitude slightly below the linear prediction ($\varOmega =1$), and with the response curve marginally asymmetric (see figure 13a). Both features are perfectly captured by the Duffing solution. For increasing driving amplitude hysteresis emerges, and the onset of hysteresis is in agreement with the analytical nonlinearity threshold for the Duffing oscillator ($A_c\approx 0.13\,\%$, see Brennan et al. Reference Brennan, Kovacic, Carrella and Waters2008). In the lower branch, the waves remain quasi-planar (as is the case for all $A$), whereas in the upper branch the wave peaks becomes higher and slightly more pointed. All our analyses suggest that the flow behaves here like an ideal Duffing oscillator, which is also reflected in the high quality of the fit (with discrepancies at the level of measurement uncertainty).

At $A=0.32\,\%$ the agreement with the fit is still excellent, but in the experiments the jump-down transition occurs much earlier than in the Duffing curve. This behaviour becomes more pronounced at $A=0.64\,\%$, whereas the agreement in the jump-up transition point ($\varOmega \approx 0.95$) remains excellent. In the Duffing oscillator, the jump-up transition depends on the nonlinearity only, whereas the jump-down transition depends also on the damping (Brennan et al. Reference Brennan, Kovacic, Carrella and Waters2008). The premature jump-down transition in our experiments suggests a nonlinear increase of the damping, which is in line with the emergence of highly dissipative waves (with liquid steeply rising at the wall or wave breaking, see figure 10c,d) close to resonance.

As shown earlier, a more useful characterisation of the deviations from Duffing dynamics is provided the phase lag between excitation and response (shown in figure 13b). For all driving amplitudes, we find that the phase lag deviates from the Duffing prediction well before the resonance is reached. At the resonance (and thus at the jump-down transition) the phase of the sloshing always lags the driving by $90^{\circ }$. Interestingly, in Duffing oscillators this jump-down transition at the response maxima occurs exactly at $90^{\circ }$ and it seems that this feature persists in our experiments despite the complex flow dynamics and the increased damping.

7.2. Comparison of the experiments to the response curves of the multimodal model

In this section, we quantitatively compare our experimentally obtained response curves to the analytical predictions of the multimodal model. We recall that the damping rate was determined from experimental measurements (see § 3.3), as for the Duffing equation, whereas all other parameters are predicted by the theory and cannot be freely fitted. As shown earlier, the horizontal motion of the liquid's centre of mass $X$ provides an unambiguous determination of the sloshing amplitude and phase from experimental measurements, and is less sensitive to the exact wave shape than the surface elevation at a particular location. The horizontal oscillation of the liquid's centre of mass $X$ can be computed from the multimodal model as

(7.1)\begin{equation} X_N(t) = \frac{1}{w\,h} \int_{{-}w/2}^{w/2} \xi\left[ \sum_{n=1}^{N} (w \beta_n(t) f_n(\xi) ) + h\right] \,{\textrm{d}} \xi /w. \end{equation}

The response curve of the first three modes ($N=3$) is shown together with the experimental measurements in figure 14.

Figure 14. Comparison between the response curves from multimodal model (lines) and the experimentally measured values (symbols). (a) Response of the centre of mass amplitude for different forcing amplitudes $A$. For the curves from the model the centre of mass location is computed using (7.1) including the first three modes $N=3$. The amplitude $\hat X$ is selected as the maximum displacement over one oscillation period $t \in [0,T]$. (b) Response curve of the phase lag between sloshing oscillation and excitation from the multimodal model with (5.11) and from the experiment.

At the lowest driving amplitude $A=0.09\,\%$ the model agrees well with the measurements sufficiently far from the resonance. However, in its vicinity no hysteresis is observed in the experiments and the response maximum reaches only approximately half the value of the model prediction. We stress that at these parameters the flow dynamics is almost linear (with quasi-planar surface waves and frequency spectra similar to the one shown exemplarily in figure 10(i)). For increasing driving amplitudes $A$ the lower branch and its jump-up transitions are well described in the multimodal model. Here, the wave shape is almost harmonic (as in figure 8a) leading to the excellent agreement of the response. At the upper branch, the sloshing amplitude obtained from the experiments is systematically below the model prediction (between $8\,\%$ and $30\,\%$) with growing deviations as the nonlinear response maxima are approached.

The damping rate is the only parameter not derived rigorously in the multimodal model. Modifying this value has little influence on the sloshing regimes, where the model predictions agree well with the experiments (e.g. at the lower branch solution and at the jump-up transitions). By contrast, the value of the damping ratio significantly affects the resonance amplitudes and the occurrence and the size of the hysteresis. If the damping ratio is decreased (e.g. $\gamma = 5.57 \times 10^{-3}$ from Keulegan's theory), the discrepancies with the experiment worsen. On the other hand, selecting larger values of the damping ratio (e.g. $\gamma = 15 \times 10^{-3}$), leads to better predictions of the hysteresis regions, especially at low amplitudes ($A=0.09\,\%, 0.16\,\%$). The prediction of high-amplitude sloshing remains, however, poor independently of the value of $\gamma$. These observations indicate that the assumption of linear damping (applied only to the first mode) is at the root of causing the deviations from the experimental observations. This is further confirmed when comparing the phase difference $\phi$ between the driving and the sloshing response, shown in figure 14(b). Overall discrepancies are very large, even in regimes in which the sloshing amplitude $\hat X$ is reasonably predicted, and increase with the driving amplitude. Further comparsion of the surface displacement amplitude from the multimodal model and the measurements is given in Appendix A.

8.1. Potential flow theory: multimodal model

The multimodal model developed by Faltinsen & Timokha (Reference Faltinsen and Timokha2009) is based on potential flow theory. It is analytical and does not have free fitting parameters. It is able to predict the full free-surface elevation and response curves for varying filling ratios, where for example the resonance shifts from softening to hardening. The multimodal analysis applied here is a single dominant system with three degrees of freedom and has therefore a solid physical foundation. Intrinsically, dissipation and damping are not included because of the underlying potential-flow equations. Typically, a linear viscous damping is added to the equation for the first mode only. The corresponding damping coefficient can either be experimentally determined or analytically estimated following Keulegan (Reference Keulegan1959). In our study, the experimentally determined coefficient leads to a better agreement of the multimodal response curves with the experimentally measured ones. In general, the model only gives accurate amplitude predictions on the lower branch or far the resonance. The sloshing amplitude at the resonance peak and the size of the hysteresis band are strongly overestimated even for small driving. The predicted phase difference between driving and response substantially deviates from the experimental measurements.

Our measurements indicate that the assumption of linear damping applied to the first mode only is a plausible cause of the observed deviations. As the driving increases, the amplitude of the second and third modes increases and the damping associated with them likely become non-negligible. We, however, stress that the multimodal model provides the best purely analytical, nonlinear prediction of the sloshing behaviour. It is, for example, able to describe the competition and interaction between modes in shallow systems or in other tank geometries (e.g. of cylindrical form). Our rectangular tank geometry and large filling height, on the other hand, were selected to generate a sloshing dynamics with a well-distinguished response peak (i.e. free of mode competition) and a quasi-two-dimensional flow. This simple setting paves the way for modelling the dynamics with the two following approaches.

8.2. Mechanical mass-spring model: Duffing oscillator

Simple mass-spring models have long been used in aerospace engineering to describe sloshing (Dodge Reference Dodge2000). The (linear) harmonic oscillator emerges from a mechanical analogy, which enables the analytical calculation of the equivalent forcing amplitude and the equivalent spring constant. The Duffing equation extends the harmonic oscillator by including a nonlinear (cubic) spring deformation and therefore lacks an a priori justification for modelling sloshing. A surprisingly good agreement between response curves of the amplitudes is achieved with a single free fitting parameter (e.g. the nonlinearity constant $\epsilon$). Measuring a single response curve in the experiments is sufficient to obtain $\epsilon$, and therefore also the response curves at other driving amplitudes (including the correct threshold value for the resonance to become hysteretic). However, at large driving amplitude the height of the resonance peak and the size of the hysteresis are overestimated. The predicted phase lag between driving and response agrees well at low sloshing amplitudes and deviates progressively for increasing amplitude. Although the agreement of the experiments with the Duffing oscillator is overall excellent, it is likely to fail for a further increasing driving amplitude, low filling height or other tank geometries. These would lead to highly complex effects and internal resonances that cannot be described by the Duffing equation anymore. Our study and the discovery of the Duffing dynamics might open avenues for modelling approaches that have the potential to capture and predict such effects. One of the most promising approaches is described in the following.

8.3. Dynamical systems approach: spectral submanifold theory

Sloshing is described by an evolutionary partial differential equation (the Navier–Stokes equation), i.e. formally by an infinite-dimensional dynamical system (Ockendon & Ockendon Reference Ockendon and Ockendon1973). The main idea of the spectral submanifold theory is that the relevant dynamics may be lower-dimensional and limited to a spectral submanifold. The excellent agreement of our experiments (at low to moderate driving amplitudes) with the forced–damped Duffing equation can only be explained by such a behaviour. It seems that the sloshing dynamics of our system is confined to a low-dimensional attractor on which the leading-order dynamics is just the Duffing equation. The recently introduced concept of time-periodic, spectral submanifolds guarantees in these types of oscillatory problems precisely a two-dimensional, time-periodic, attracting invariant manifold (Breunung & Haller Reference Breunung and Haller2018). The increasing deviations from the Duffing equation at larger driving amplitudes indicate that higher-order polynomial models are required to approximate the dynamics on the reduced two-dimensional, time-periodic spectral submanifold. In principle these polynomials could be fitted to the experimental data, but this approach is outside of the scope of this work. Note that in our experiment internal resonances seem to be absent, although an exception could be the reported period-three motion. These are likely to occur in other tank geometries or at lower filling levels and would lead to a dynamics that cannot be described anymore in a two-dimensional spectral submanifold. In that case, higher-dimensional spectral submanifolds (with coupled Duffing-type equations) might provide again accurate predictions.