2. Experimental set-up

The tank used for the experimental campaign at CNR-INM is $L=1$ m long and tall, $D=0.1$ m wide and filled with water up to $h_{FS}=0.35$ m. The total mass of the liquid is therefore equal to $m_l=34.76$ Kg. Based on the filling height $h_{FS}$ selected, the natural sloshing periods can be derived from (Faltinsen & Timokha Reference Faltinsen and Timokha2009)

(2.1)\begin{equation} T_n = \displaystyle{\frac{2 {\rm \pi}}{\sqrt{\displaystyle{\dfrac{g n {\rm \pi}}{L}} \tanh{\displaystyle{\dfrac{n {\rm \pi}h_{FS}}{L}}}}}} \quad n=1, 2, \ldots \end{equation}

by setting $n$ equal to the selected mode. For example, by considering that $g$ is the gravitational acceleration, the period of the first mode is $T_1=1.265$ s.

To guarantee a purely sinusoidal lateral movement, $x(t)=A \sin {(2 {\rm \pi}t/T)}$, an ad hoc mechanical system was designed; $A$ and $T$ are the amplitude and period of the prescribed motion, respectively. The small breadth of the tank, i.e. $D=0.1$ m, ensures an almost 2-D flow in the sloshing plane. Five capacitive wire probes are positioned inside the tank. The first two are located at a distance of $1$ cm and $5$ cm from the left side. Two other probes are positioned symmetrically on the right side, while the last wire probe is located in the centre of the tank (i.e. $50$ cm from both sides of the tank). The corresponding dimensionless measurements of the surface elevations $h_i$ are indicated as $\eta _i = (h_i-h_{FS})/L$ with $i=1,5,50,95,99$.

During the tests, flow visualizations are obtained through a digital video camera JAI CV-M2. This camera has a spatial resolution of $1600 \times 1200$ pixels and a frame frequency of 15 Hz. It is positioned in front of the tank, far enough away from it to record the entire flow pattern. A wire potentiometer was used to evaluate the position of the tank. A synchronizer is used to trigger the start of all acquisition systems, which are characterized by different sampling rates. A sketch of the experimental set-up is given in figure 1.

Figure 1. Experimental box sketch with highlight of the filled volume. (a) A photo of the experimental set-up; (b) a sketch of the experimental arrangement of the probes.

Multiple tests were performed with the same parameters to verify the repeatability of the results. The wave heights measured by the capacitive wire probes were also verified with those extracted from the videos in the time windows in which the digital images were archived. Since this is a medium-sized experimental apparatus, no critical issues emerged in the movement regimes studied. The mechanical limits are mainly related to the range of oscillation frequencies allowed by the electric motor and mechanical guide used to enforce the movement of the tank.

The results of the CNR-INM experimental campaign were also published in the book by Faltinsen & Timokha (Reference Faltinsen and Timokha2021). In the present work, the experimental data were used in synergy with the numerical outcomes for a detailed analysis of peculiar nonlinear sloshing regimes.

3. Numerical solver

In the present work a 2-D fluid domain $\varOmega$ is considered with boundaries which are composed of the tank walls $\partial \varOmega _B$ and the free surface $\partial \varOmega _F$. Only the liquid phase is considered and modelled as a barotropic weakly compressible medium. The tank moves along the $x$-axis and the equations are formulated in the non-inertial frame of reference (Ni-FoR). According to these assumptions, the flow evolution is described as

(3.1) \begin{equation} \left. \begin{array}{@{}ll@{}} \displaystyle \dfrac{{\rm D} \rho}{{\rm D} t} ={-}\rho \,\mbox{div}( \boldsymbol{u}), & \displaystyle \dfrac{D \boldsymbol{u}}{D t} = \dfrac{\mbox{div} (\mathbb{T})}{\rho} +\boldsymbol{g}-a_{tank}(t)\,\boldsymbol{e}_1\\ \displaystyle \dfrac{{\rm D} e}{{\rm D} t} = \dfrac{\mathbb{T}:\mathbb{D}}{\rho}, & \displaystyle \dfrac{{\rm D} \boldsymbol{r}}{{\rm D}t} = \boldsymbol{u}, \quad p = f(\rho), \end{array} \right\} \end{equation}

with ${\rm D}/{\rm D}t$ the Lagrangian derivative, $\boldsymbol {r}$ the local position, $\boldsymbol {u}$ the flow velocity, $\rho$ the fluid density, $e$ the specific internal energy, $\mathbb {T}$ the stress tensor, $\mathbb {D}$ the rate of strain tensor, $\boldsymbol {g}$ the acceleration related to the gravitational field, $a_{tank}(t)$ the tank acceleration and $\boldsymbol {e}_1$ the unit vector of the $x$-axis,

The fluid is assumed Newtonian and the flow isothermal, while the surface tension is neglected, i.e. $\mathbb {T}=[- p+\lambda \,\mbox {div}(\boldsymbol {u})]\,\mathbb {I}+2 \mu \mathbb {D}$, where $\mu$ and $\lambda$ are the primary and secondary dynamic viscosities of the liquid and $\mathbb {I}$ the identity tensor.

Within the single-phase model, all the physics related to the air entrapment in the fluid domain are neglected. This assumption may appear inappropriate for violent sloshing simulations. However, in Marrone et al. (Reference Marrone, Colagrossi, Di Mascio and Le Touzé2016) and Malan et al. (Reference Malan, Pilloton, Colagrossi and Malan2022), it was shown that the energy dissipation in violent flows can be evaluated with sufficient precision even in the single-phase hypothesis. Further confirmation can also be found in Bouscasse et al. (Reference Bouscasse, Colagrossi, Souto-Iglesias and Cercos-Pita2014a,Reference Bouscasse, Colagrossi, Souto-Iglesias and Cercos-Pitab), where the study of mechanical energy dissipation induced by sloshing and wave breaking in a fully coupled angular motion system was investigated. In this work, the single-phase model was shown to be able to predict experimental results with reasonable accuracy.

By considering that the temperature is assumed constant, since the effects of its variation are negligible with good approximation, it is assumed that the pressure $p$ depends on the density, only. Furthermore, since a weakly compressible condition is assumed, a simple linear equation of state can be considered

(3.2)\begin{equation} p = c_0^2 ( \rho - \rho_{0} ), \end{equation}

where $c_0$ plays the role of a constant speed of sound of the liquid and $\rho _0$ is the density at the free surface (where $p$ is assumed to be equal to zero).

The weakly compressible hypothesis implies the following constraint:

(3.3)\begin{equation} c_0 \gg \max{\left(\Delta U_{max},\sqrt{\frac{\Delta p_{max}}{\rho_0}}\right)}, \end{equation}

where $\Delta U_{max}$ and $\Delta p_{max}$ are, respectively, the maximum fluid speed and the maximum pressure variations expected in $\varOmega$ during the time evolution. Considering that the temporal integration is carried out with a time step linked to the value of $c_0$, the latter is always set lower than its physical counterpart. Constraint (3.3), however, must be verified to ensure compliance with the weakly compressible regime.

In the present paper, the whole set of numerical simulations spans from $Re=3\times 10^4$ to $Re=3 \times 10^5$, where the Reynolds number is defined as

(3.4)\begin{equation} Re = \displaystyle{\frac{2 {\rm \pi}A}{T}} \displaystyle{\frac{L}{\nu_w}}, \end{equation}

where $\nu _w = 10^{-6} \ {\rm m}^2\ {\rm s}^{-1}$ is the kinematic viscosity of the water, $A$ is the amplitude of tank oscillation, $T$ the excitation period and $L$ the tank length (1 m in the present investigations). By considering the Reynolds numbers involved, a sub-grid model is needed. Similarly to the articles of Christensen & Deigaard (Reference Christensen and Deigaard2001) and Christensen (Reference Christensen2006), a LES with a classic Smagorinsky model is adopted for simulating breaking waves.

The LES model is adopted in the SPH framework by several authors see e.g. Lo & Shao (Reference Lo and Shao2002), Rogers & Dalrymple (Reference Rogers and Dalrymple2005), Violeau & Issa (Reference Violeau and Issa2007), Di Mascio et al. (Reference Di Mascio, Antuono, Colagrossi and Marrone2017) and Meringolo et al. (Reference Meringolo, Marrone, Colagrossi and Liu2019). In this work, the $\delta$-LES-SPH model of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021a) is used, the peculiarity of which is that the LES model was introduced within a quasi-Lagrangian formalism.

3.1. Brief recall of the $\delta$-LES-SPH model

Within the SPH model the fluid domain $\varOmega$ is discretized in a finite number of particles. The differential equations for the motion of those fluid particles come from the discretization of the governing equations (3.1) according to the $\delta$-LES-SPH model

(3.5) \begin{equation} \left. \begin{array}{@{}l@{}} \displaystyle{\dfrac{{\rm d}\rho_i}{{\rm d} t}} = \displaystyle{\sum_j} \left[ -\rho_i (\boldsymbol{u}_{ji} + \delta \boldsymbol{u}_{ji}) + (\rho_j \delta \boldsymbol{u}_j + \rho_i \delta \boldsymbol{u}_i) \right]\boldsymbol{\cdot} \boldsymbol{\nabla}_i W_{ij} V_j + \mathcal{D}^{\rho}_i\\ \displaystyle{\dfrac{{\rm d} \boldsymbol{u}_i}{{\rm d} t}} =\displaystyle \dfrac{1}{\rho_i}\,{\sum_j} \left[-(p_j+p_i) \,\mathbb{I} + \rho_0 (\boldsymbol{u}_j \otimes \delta \boldsymbol{u}_j + \boldsymbol{u}_i \otimes \delta \boldsymbol{u}_i) \right] \boldsymbol{\cdot} \boldsymbol{\nabla}_i W_{ij} V_j\\ \quad\qquad +\, \boldsymbol{F}^v_i + \boldsymbol{g} -a_{tank}(t)\,\boldsymbol{e}_1\\ \displaystyle{\dfrac{{\rm d} \boldsymbol{r}_i}{{\rm d} t}} = \boldsymbol{u}_i+ \delta \boldsymbol{u}_i, \quad V_i(t)=m_i/\rho_i(t),\quad p=c_0^2(\rho-\rho_0), \end{array} \right\} \end{equation}

where the index $i$ refers to the considered particle and $j$ refers to neighbouring particles of $i$. The vector $\boldsymbol {F}^v_i$ is the net viscous force acting on the particle $i$. The notation $\boldsymbol {u}_{ji}$ in (3.5) indicates the differences $(\boldsymbol {u}_j-\boldsymbol {u}_i)$ and the same holds for $\delta \boldsymbol {u}_{ji}$ and $\boldsymbol {r}_{ji}$. The spatial gradients are approximated through the convolution with a kernel function $W_{ij}$. Following Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021a), a C2-Wendland kernel is adopted in the present work.

In order to recover a regular spatial distribution of particles and consequently an accurate approximation of the SPH operators (Quinlan, Lastiwka & Basa Reference Quinlan, Lastiwka and Basa2006; Nestor et al. Reference Nestor, Basa, Lastiwka and Quinlan2009), a particle shifting technique is used (see also e.g. Lind et al. Reference Lind, Xu, Stansby and Rogers2012). For the sake of brevity, the specific expression of the shifting velocity $\delta \boldsymbol {u}$ is not reported here, but the interested reader may refer to Marrone et al. (Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez2021b) and Michel et al. (Reference Michel, Durante, Colagrossi and Marrone2022), where violent sloshing problems were investigated.

The time derivative ${\rm d}/{\rm d}t$ used in (3.5) indicates a quasi-Lagrangian derivative, i.e.

(3.6)\begin{equation} \frac{{\rm d}({\bullet})}{{\rm d}t}:=\frac{\partial ({\bullet})}{\partial t}+\boldsymbol{\nabla}({\bullet}) \boldsymbol{\cdot} (\boldsymbol{u}+\delta \boldsymbol{u}), \end{equation}

since the particles are moving with the modified velocity $(\boldsymbol {u}+\delta \boldsymbol {u})$ and the first two equations of (3.5) are written following an arbitrary Lagrangian–Eulerian approach. Due to this, the continuity and the momentum equations contain terms with spatial derivatives of $\delta \boldsymbol {u}$ (for details, the interested reader is referred to Antuono et al. Reference Antuono, Sun, Marrone and Colagrossi2021b).

The mass $m_i$ of the $i$th particle is assumed to be constant during its motion. The particles are set initially on a Cartesian lattice with spacing $\Delta r$, and hence, the volumes $V_i$ are initially set as $\Delta r^2$. The particle masses $m_i$ are calculated through the initial density field (using the equation of state and the initial pressure field). While the particle masses $m_i$ remain constant during the time evolution, the volumes $V_i$ change over time in accordance with the particle density (see bottom line of (3.5)).

To avoid instability in the pressure field, the diffusive term $\mathcal {D}^{\rho }_i$, introduced by Antuono, Colagrossi & Marrone (Reference Antuono, Colagrossi and Marrone2012b), is added into the continuity equation. For brevity, $\mathcal {D}^{\rho }_i$ is not reported here, the interested reader can also find more details in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021a) and in Meringolo et al. (Reference Meringolo, Marrone, Colagrossi and Liu2019), where the intensity of this term is determined dynamically in space and time.

The viscous forces $\boldsymbol {F}^v$ are expressed as

(3.7) \begin{equation} \left.\begin{array}{@{}ll@{}} \boldsymbol{F}^v_i:=\displaystyle K\,\sum_j(\mu+\mu^T_{ij}){\rm \pi}_{ij}\boldsymbol{\nabla}_i W_{ij} V_j & K:=2(n+2)\\ {\rm \pi}_{ij} :=\displaystyle \dfrac{\boldsymbol{u}_{ij}\boldsymbol{\cdot}\boldsymbol{r}_{ij}}{||\boldsymbol{r}_{ji}||^2} \quad \mu_{ij}^T:=2\dfrac{\mu_i^T\mu_j^T}{\mu_i^T+\mu_j^T} & \mu_i^T:=\rho_0(C_S\,l)^2\,||\mathbb{D}_i||, \end{array}\right\} \end{equation}

where $n$ is the number of spatial dimensions, $l=4 \Delta r$ is the radius of the support of the kernel $W$ for two spatial dimensions and represents the length scale of the filter adopted for the LES sub-grid model; $C_S$ is the so-called Smagorinsky constant, set equal to 0.18 (as in Smagorinsky Reference Smagorinsky1963; Bailly & Comte-Bellot Reference Bailly and Comte-Bellot2015) and $||\mathbb {D}||$ is a rescaled Frobenius norm, namely $||\mathbb {D}||=\sqrt {2\mathbb {D}:\mathbb {D}}$. The viscous term (3.7) contains both the effect of the physical viscosity $\mu$ as well as of the turbulent stresses $\mu _i^T$. In order to dump the turbulent eddies near the wall boundaries, a classical van Driest damping function is employed (Van Driest Reference Van Driest1956) (for more details, see also Pilloton et al. Reference Pilloton, Michel, Colagrossi and Marrone2023).

A fourth-order Runge–Kutta scheme is adopted to integrate in time the system (3.5). The time step $\Delta t$ is obtained as the lowest among the following constraints, related to the Courant–Friedrichs–Lewy conditions

(3.8) \begin{equation} \left.\begin{array}{@{}lll@{}} \displaystyle\Delta t_v = 0.5 \min_i \dfrac{\Delta r^2\rho_i}{(\mu+\mu_i^T)} , & \Delta t_a = 0.3 \min_i \sqrt{ \dfrac{\Delta r}{\| \boldsymbol{a}_i \|}} , & \Delta t_c = 2.0 \left(\dfrac{\Delta r}{c_0}\right)\\ \Delta t=\min(\Delta t_v, \Delta t_a, \Delta t_c), & & \end{array}\right\} \end{equation}

where $\| \boldsymbol {a}_i \|$ is the particle acceleration, $\Delta t_v$ is the time step related to viscosity, $\Delta t_a$ is the advective time step and $\Delta t_c$ is the acoustic time step (see e.g. Colagrossi et al. Reference Colagrossi, Rossi, Marrone and Le Touzé2016). For the cases studied in this work, the last two constraints (involving $\Delta t_a$ and $\Delta t_c$) are always the most critical.

3.3. Evaluation of the slosh dissipation

Following the analysis performed in Marrone et al. (Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez2021b) and Michel et al. (Reference Michel, Durante, Colagrossi and Marrone2022), the $\delta$-LES-SPH energy balance, in terms of power, can be written as

(3.9) \begin{equation} \left.\begin{array}{@{}l@{}} \displaystyle\dot{\mathcal{E}}_K+\dot{\mathcal{E}}_P -\mathcal{P}_{NF} = \mathcal{P}_V + \mathcal{P}_V^{turb} + \mathcal{P}_N \\ \mathcal{E}_K(t)=\dfrac{1}{2}\sum\limits_im_iu_i^2,\quad \mathcal{E}_P(t)=\sum\limits_im_i gy_i ,\quad \mathcal{P}_{NF}=\sum\limits_im_ia_{tank}(t)\,\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{u}_i, \end{array}\right\} \end{equation}

where, on the left-hand side, $\mathcal {E}_K$ and $\mathcal {E}_P$ are the kinetic and potential energy of the particle system in the Ni-FoR moving with the tank. The vertical position of the generic $i$th particle is indicated with $y_i$. Finally, $\mathcal {P}_{NF}$ is the power related to non-inertial forces. The potential energy related to the compressibility of the liquid is negligible under the weakly compressible assumption, therefore it is not considered in the energy balance (for more details see Antuono et al. Reference Antuono, Marrone, Colagrossi and Bouscasse2015).

The right-hand side of the energy balance (3.9) contains the dissipation term related to the physical viscosity $\mathcal {P}_V$ and to the eddy viscosity $\mathcal {P}^{turb}_V$, while $\mathcal {P}_N$ takes into account the numerical effect of the density diffusion and of the particle shifting $\delta \boldsymbol {u}$ (see also Michel et al. Reference Michel, Antuono, Oger and Marrone2023; Sun et al. Reference Sun, Pilloton, Antuono and Colagrossi2023). The power related to the viscous forces is directly evaluated through the expression (3.7)

(3.10)\begin{equation} \mathcal{P}_V + \mathcal{P}_V^{turb} = 5 \sum_i\sum_j (\mu + \mu_{ij}^T) {\rm \pi}_{ij} \boldsymbol{u}_{ij} \boldsymbol{\cdot} \boldsymbol{\nabla}_i W_{ij} V_i V_j, \end{equation}

where the quantity $\mathcal {P}_V^{turb}$ refers to the viscous dissipation of the modelled sub-grid scales.

The energy dissipated by the fluid is then evaluated by integrating in time equation (3.9)

(3.11a)\begin{equation} [\mathcal{E}_K+\mathcal{E}_P](t) - [\mathcal{E}_K+\mathcal{E}_P](t_0) - \mathcal{W}_{NF}(t) = \mathcal{E}_{diss}(t)\end{equation}

(3.11b)\begin{equation} \mathcal{W}_{NF}(t) = \displaystyle\int\nolimits_{t_0}^t \left[\sum\limits_i m_i ({-}a_{tank}(t)\,\boldsymbol{e}_1 \boldsymbol{\cdot}\boldsymbol{u}_i) \right] {\rm d}t, \quad \mathcal{E}_{diss}(t) = \displaystyle\int\nolimits_{t_0}^t \left(\mathcal{P}_V + \mathcal{P}_V^{turb} + \mathcal{P}_N \right) {\rm d}t, \end{equation}

where $\mathcal {W}_{NF}$ is the work performed by the non-inertial forces and $[\mathcal {E}_K+\mathcal {E}_P](t_0)$ is the mechanical energy at the initial instant $t_0$.

The energy dissipated by the fluid can be either evaluated by the left-hand side of top equation (3.11a,b), or by the direct estimate of $\mathcal {E}_{diss}$ with the bottom expression in (3.11a,b). As performed in Malan et al. (Reference Malan, Pilloton, Colagrossi and Malan2022) and in Marrone et al. (Reference Marrone, Saltari, Michel and Mastroddi2023), both these approaches can be adopted in order to verify that a specific SPH model is able to close the energy balance accurately.

Considering that the fluid is initially at rest, the energy balance (3.11a,b) can be reshaped as

(3.12) \begin{equation} \left.\begin{array}{@{}l} \displaystyle\mathcal{E}_K(t)+m_lg\left[y_G(t)-y_G(t_0)\right]+m_l\int_{t_0}^t\dot{x}_G(t)\,a_{tank}(t)\,{\rm d}t =\mathcal{E}_{diss}(t)\\ \displaystyle x_G(t):=\dfrac{\sum_im_ix_i(t)}{m_l} \quad y_G(t):=\dfrac{\sum_im_iy_i(t)}{m_l}, \end{array} \right\} \end{equation}

where $x_G$ and $y_G$ are the horizontal and vertical coordinates of the liquid centre of mass in the Ni-FoR, with $m_l=\sum _i m_i$ the total mass of the liquid inside the tank.

Assuming that the tank oscillates with a harmonic law, after a sufficiently long transient a periodic or quasi-periodic regime can be attained. It may take several periods to reach this condition, due to the highly nonlinear behaviour of the sloshing phenomenon. In the ‘quasi-periodic’ condition, the mechanical energy becomes negligible compared with the work $\mathcal {W}_{NF}$, which almost coincides with the dissipated energy.

In particular, the energy dissipated during the $k$th period can be evaluated as

(3.13)\begin{equation} \mathcal{E}_{diss}^{(k)} = m_l\int_{kT}^{kT+T} \dot{x}_G(\tau)\,a_{tank}(\tau) \,{\rm d}\tau. \end{equation}

As shown in § 8, after the transitory stage lasting $N_{P0}$ periods, $\mathcal {E}_{diss}^{(k)}$ remains almost constant for $k>N_{P0}$, hence, it is possible to associate an average dissipation power during $N_p$ periods as

(3.14)\begin{equation} \bar{\mathcal{P}}_{diss}=\frac{m_l}{N_PT} \int_{N_{P0}T}^{(N_{P0}+N_P)T} \dot{x}_G(\tau)\,a_{tank}(\tau) \,{\rm d}\tau. \end{equation}

The above equation highlights the role of the motion of the centre of mass of the fluid in the Ni-FoR concerning the energy dissipated by the liquid sloshing. The phase lag between the horizontal tank motion $x_{tank}(t)$ and the fluid centre of mass $x_G(t)$ plays a crucial role in the slosh dissipation (see e.g. Marrone et al. Reference Marrone, Colagrossi, Calderon-Sanchez and Martinez-Carrascal2021a,Reference Marrone, Colagrossi, Gambioli and González-Gutiérrezb; Malan et al. Reference Malan, Pilloton, Colagrossi and Malan2022). As is clear from (3.14), high dissipation is obtained when $\dot {x}_G(t)$ and $a_{tank}(t)$ are in anti-phase, whereas low dissipation occurs when they are in quadrature (see also Malan et al. Reference Malan, Pilloton, Colagrossi and Malan2022; Saltari et al. Reference Saltari, Pizzoli, Coppotelli, Gambioli, Cooper and Mastroddi2022; Marrone et al. Reference Marrone, Saltari, Michel and Mastroddi2023).

5. Tripling-period and doubling-frequency modes analysed with the Hilbert–Huang transform

In the present section the Hilbert–Huang transform (HHT) is used to make some aspects related to doubling and tripling modes more clear. The HHT was designed for nonlinear and non-stationary signals (Huang et al. Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998) and it works by two successive steps: the empirical mode decomposition (EMD) and the Hilbert transform (HT). The EMD is used to decompose a signal into a number of intrinsic mode functions (IMFs). An IMF is an oscillatory function with time-varying frequencies that represents the characteristics of non-stationary signals (Lee et al. Reference Lee, Chang, Hsieh, Deng and Sun2012).

In the time range $[t_1-t_2]$ the reference signal $X(t)$ is expressed as a combination of IMFs as

(5.1)\begin{equation} X(t) = \displaystyle{\sum_{i=1}^n} \,\mbox{IMF}_i(t) + \mathcal{R}_n(t), \end{equation}

where $n$ is the total number of IMFs and $\mathcal {R}_n(t)$ is the residue function. It is worth noting that, from a theoretical point of view, IMFs are not strictly orthogonal, although this property is actually satisfied, as verified a posteriori.

The HT $\mathcal {H}$ of IMF$(t)$ consists of a convolution with the Cauchy kernel $k(t)=1/{\rm \pi} t$

(5.2)\begin{equation} \mathcal{H}(\mbox{IMF}_k(t))={\frac {1}{{\rm \pi} }} \,{\int\hskip -1,05em -\,} _{-\infty }^{+\infty }{\frac {\mbox{IMF}_k(\tau )}{t-\tau }}\,\mathrm {d} \tau \quad\mbox{with} \ k = 1, \ldots, n, \end{equation}

where the principal value of the integral is considered. The signal $X(t)$ is represented in the complex plane with the function $\boldsymbol {Z}(t)$ (complex quantities are indicated in bold) defined as

(5.3)\begin{equation} \left.\begin{gathered} \boldsymbol{Z}_k(t) = \mbox{IMF}_k(t) + \boldsymbol{i} \mathcal{H}(\mbox{IMF}_k(t))= A_k(t) \exp{[ \boldsymbol{i} \theta_k(t) ]} \quad \mbox{with} \ k = 1, \ldots, n \\ A(t) = |\boldsymbol{Z}(t)| , \quad f(t) = \displaystyle{\frac{1}{2 {\rm \pi}} \,\frac{{\rm d}\theta(t)}{{\rm d}t}} \quad \mbox{with} \ \theta(t) = \angle{\boldsymbol{Z}(t)} \end{gathered}\right\}, \end{equation}

with $\boldsymbol {i} = \sqrt {-1}$, $A(t)$ the instantaneous amplitude and $f(t)$ the instantaneous frequency of $\boldsymbol {Z}(t)$.

The time behaviours of the frequencies of different IMFs, made non-dimensional with the excitation frequency $f_e=1/T$, are drawn in figure 15 for $T=1.304 T_1$ and in figure 17 for $T=0.867 T_1$. These cases correspond to a doubling-frequency and to a tripling-period bifurcations, respectively. In order to highlight the most energetic modes, the curves are contoured with the non-dimensional amplitude $A(t)/L$. The ordinate $f^*=fT$ refers to dimensionless frequencies, so the dimensionless excitation frequency $f_eT$ is 1 by definition. Considering that the numbering of IMF modes is arbitrary, we decided to number the modes by their energy content, so that the mode IMF$_1$ is the most energetic and so on.

Figure 15. Hilbert transform for the case $A=0.03 L$ and $T=1.304 T_1$. The different lines correspond to the time trends of the frequencies related to the different modes. The lines are contoured with the amplitudes of the corresponding modes for varying $t$. The frequencies are non-dimensional with the exciting carrier $1/T$.

Since the energy of $X(t)$ is in general roughly equal to the sum of the energies of its modes, the distribution of the signal energy among IMF modes is discussed. It is worth underlining that, with a little abuse of language, we will refer in the following to the mode energy, although we are actually dealing with its amplitude.

Figure 15 represents the time trends of IMF mode frequencies at $T=1.304 T_1$, where a doubling-frequency scenario occurs. Although the system is excited at frequency $f_e$, a significant part of the energy moves to $2 f_e$, making this mode comparable, in terms of energy, to the first harmonic. This happens because the second harmonic of the excitation frequency is close to the second resonance, i.e. $2f_e \sim f_2=1/T_2$, as indicated by the formula (2.1). This induces an energy transfer toward $2 f_e$, similarly to what observed by Wu (Reference Wu2007), who performed similar investigations with a linear perturbation theory. The second mode is excited by nonlinearity so that, in figure 15, it is possible to appreciate that the mode IMF$_1$, related to frequency $2f_e$, starts at approximately 10 oscillation periods and it persists over time without appreciable perturbations. With analogous arguments, modes IMF$_3$ and IMF$_4$, linked to frequencies $3 f_e$ and $4 f_e$, become rather similar as well. The final emerging pattern is typical of a doubling-frequency regime. Figure 16 shows the time histories of the sum of the first two IMFs ${\rm IMF}_1+{\rm IMF}_2$. The plot depicts the onset and the persistence of the doubling-frequency regime.

Figure 16. Sloshing case $A=0.03 L$ and $T=1.304 T_1$: time histories of the sum of the first two IMFs ${\rm IMF}_1+{\rm IMF}_2$ calculated from the $\delta$-LES-SPH $\eta _5$ time signal.

The case $T= 0.867 T_1$, where a tripling-period bifurcation occurs, is more complex. As shown in figure 17, the mode ${\rm IMF}_1$ related to the excitation frequency $f_e$ rapidly appears and persists as almost constant in time. The mode ${\rm IMF}_4$, related to first resonance $f_1 = 1/T_1$, arises during the early periods of oscillation as an effect of the initial ramp used for the simulation but, after roughly $40$ oscillations, the mode is almost totally damped with an amplitude that becomes negligible.

Figure 17. Hilbert transform for the case $A=0.03 L$ and $T=0.867 T_1$. The different lines correspond to the time trends of the frequencies related to the different modes related to the $\delta$-LES-SPH $\eta _5$ signal. The lines are contoured with the amplitudes of the corresponding modes for varying $t$. The frequencies are made non-dimensional with the exciting carrier $1/T$.

The second and third sloshing natural frequencies (i.e. $f_2 = 1/T_2$ and $f_3 = 1/T_3$, respectively) are, in this case, related to $f_e$ as $f_2 \sim f_e + 1/3 f_e$ and $f_3 \sim f_e + 2/3 f_e$. This implies that, when the energy of ${\rm IMF}_4$ is redistributed among other frequencies, it falls on modes ${\rm IMF}_2$ and ${\rm IMF}_5$, corresponding to $f_2$ and $f_3$. These latter become energetic enough to persist in time and to become sub-harmonics between $f_e$ and $2 f_e$. The modes ${\rm IMF}_2$ and ${\rm IMF}_5$ interact with ${\rm IMF}_3$, making the modes related to sub-harmonics $2 f_e - f_2 \sim 2/3f_e$ and $2 f_e - f_3 \sim f_e/3$ to appear. The final frequency spectrum attained by the system assumes the typical peaked shape of a tripling-period one (see figure 10c).

Figure 18 shows the time histories of the sum of the first two IMFs ${\rm IMF}_1+{\rm IMF}_2$. The plot depicts the onset and the persistence of the tripling-period regime.

Figure 18. Sloshing case $A=0.03 L$ and $T=0.867 T_1$: time histories of the sum of the first two IMFs ${\rm IMF}_1+{\rm IMF}_2$ calculated from the $\delta$-LES-SPH $\eta _5$ time signal.

Finally, the case at $T= 1.022 T_1$ is depicted in figure 19, where the experimental outcomes and the numerical solution are compared. As shown in panel (a) of the figure, the analysis of the experimental result shows a dynamics very similar to $T=0.867 T_1$, discussed above. The mode corresponding to $f_1$ appears during initial oscillation periods, caused by the ramp, and rapidly damps out. The energy is then redistributed among other frequencies and the tripling-period pattern forms again, similarly to that just observed in figure 17. Conversely, the 2-D restriction assumed in the numerical simulations leads here to a discrepancy. The energy of mode related to frequency $f_1$ remains confined to its initial state and does not spread among other frequencies, although it starts to oscillate chaotically.

Figure 19. Hilbert transform for the case $A=0.03 L$ and $T=1.022 T_1$. The different lines correspond to the time trends of the frequencies related to the different modes. The lines are contoured with the amplitudes of the corresponding modes for varying $t$. The frequencies are non-dimensional with the exciting one $1/T$. (a) Experimental data, (b) $\delta$-LES-SPH outcome.

This plays the role of a chaotic modulation of the carrier signal, giving rise to a quasi-periodic pattern. For this reason, a tripling-period regime is not attained. This discordance between 2-D SPH simulations and experimental data is further discussed in § 6.

6. Tripling-period regime at $A=0.03L$ and $T= 1.022 T_1$: 2-D vs 3-D comparison

In § 5, the doubling and tripling bifurcations appearing in the WEFDs are discussed and investigated through the HHT. For $A=0.03L$ and $T = 1.022 T_1$, a discrepancy between numerical and experimental evidence is found and a tripling-period behaviour of the time signal appears from experiments, whereas numerical simulations provide a quasi-periodic mode prediction. In order to find a reason for this discrepancy, we first try to change the solid boundary condition on the tank wall. However, in this particular test case the 2-D solutions appear to be not so influenced by those changes even when a free-slip condition is considered.

The disagreement between numerical and experimental data is presumably due to the 2-D restriction, however, in order to gain a deeper insight into this problem, 3-D simulations have been carried out with a $\delta$-LES-SPH model and the results are shown in figure 20 in terms of the surface elevation time signal, and in figure 21 in terms of the free-surface evolution. The resolution used for the simulation is $N=H/\Delta r=50$, which means approximately 100 000 particles in the whole numerical domain. A free-slip boundary condition was used on the tank walls, thus avoiding the discretization of the boundary layers, in particular at low spatial resolutions (i.e. $N=50$). After many tests, this assumption was found accurate enough in two dimensions, thus it was adopted also in three dimensions.

Figure 20. Time histories of the surface elevations $\eta _5$ for the case $A=0.03 L$ and $T/T_1=1.022$. Solid line: experimental data; dashed line: 3-D $\delta$-LES-SPH simulation. Letters A,B,C and D identify the time instants for the particle displacements sketched in figure 21.

Figure 21. Free-surface evolution for the case $A=0.03 L$ and $T=1.022 T_1$ obtained with 3-D $\delta$-LES-SPH simulations.

The significant increase in computational costs did not allow a long simulation time, thus the simulation is limited to $t=100 T$ and, for the same reasons, the present investigation, including approximately 230 simulations, cannot be carried out within a 3-D environment.

Figure 20 shows a comparison between $\delta$-LES-SPH and the experimental time history of the surface elevation recorded by the probe $\eta _5$. From this plot it is possible to appreciate that, unlike the 2-D simulation, the 3-D numerical simulation is able to identify the tripling-period regime. In table 1, the comparisons among the values of $\bar {\eta }_{5 max}$ from 2-D and 3-D simulations and experiments is detailed. The average maximum of the 3-D numerical signal differs by approximately $2.5\,\%$ from the experiments, whereas the difference grows to $35.7\,\%$ with the 2-D simulation.

Table 1. Comparisons in terms of $\bar {\eta }_{5 max}$ from 2-D and 3-D simulations and from experimental outcomes. The standard deviations are also indicated in the rightmost column.

Figure 21 displays the lateral view of the particle configuration at four time instants. Particles are coloured with pressure. From this figure, it is possible to recognize the occurrence of the breaking waves which characterize the tripling-period regime, as already discussed in § 4.2.

Figure 22 shows the images recorded by the digital camera at the same instants as figure 21. The main flow features recorded by the digital camera are well predicted by the numerical solvers. As already mentioned, this pattern does not appear in two dimensions, indicating that, for this specific case, the transverse velocity component plays a key role because of the fragmentation of the free surface.

Figure 22. Digital images for the case $A=0.03 L$ and $T=1.022 T_1$ at the same time instants as figure 21.

7. Asymmetric regimes

For the problem of sloshing in a tank, Olsen (Reference Olsen1970) tacitly assumed a global symmetry between the left and right sides of the tank. Intuitively, this seems reasonable because the tank is a cuboid moving laterally leftwards and rightwards, so there is no reason for which any asymmetry should appear. For an oscillation period $T$, the right probe should record the same signal as the left one, but delayed by $T/2$. Unexpectedly, Colagrossi (Reference Colagrossi2005), Colagrossi et al. (Reference Colagrossi, Palladino, Greco, Lugni and Faltinsen2006) and Faltinsen & Timokha (Reference Faltinsen and Timokha2009) found that, at high motion amplitudes, the left $\eta _L$ and the right $\eta _R$ probes do not measure similar signals. In the present section, the presence of signal asymmetries is numerically addressed.

In order to have a quantitative evaluation of the asymmetry, a measure of the left and right probe cross-correlation is computed for each case. The correlation function is defined as (Stoica & Moses Reference Stoica and Moses2005)

(7.1)\begin{equation} R_{[ \eta_L \ \eta_R ]}(m) = \displaystyle{\sum_{k=1}^{N-m-1}} \eta_L(t_k + m \Delta t) \eta_R(t_k) \quad \mbox{with} \ t_k = t_0 + k \Delta t, \end{equation}

where $\Delta t$ is the time discretization of the signals and $m \Delta t$ is the time delay between signals. Varying the parameter $m$, the maximum value $R_{[ \eta _5 \ \eta _{95} ]}^*$ is considered. In order to get a meaningful value, the normalized cross-correlation is used

(7.2)\begin{equation} \hat{R}_{[ \eta_L \ \eta_R ]} = 1 - \displaystyle{\frac{ R_{[ \eta_L \ \eta_R ]}^* }{ \sqrt{ R_{[ \eta_L \eta_L ]}(0) \ R_{[ \eta_{R} \ \eta_{R} ]}(0) } }} . \end{equation}

For two identical signals, $\hat {R}_{[ \eta _L \ \eta _R ]}$ is equal to 0, while it tends to 1 for two non-correlated signals. As reported in figure 23 in the range $T\in (0.5 T_1, 0.7 T_1)$, the left and right probes are less correlated than in the rest of the range. In particular, for an oscillation amplitude of $A=0.03L$, the maximum asymmetry is attained for chaotic modes and specifically for the case $T=0.593 T_1$. Conversely, for $A = 0.01L$ the maximum of asymmetry is attained for quasi-periodic regimes.

Figure 23. Cross-correlation distribution with frequency for $A=0.01 L$ and $A=0.03 L$.

An example of a non-symmetric mode for $A = 0.03L$ is shown in figure 24, where the acquisitions of the free-surface elevation $\eta _1$ for the left side and $\eta _{99}$ for the right side of the tank are compared. The full signals are shown in the top panel of the figure, whereas a magnification of the range $t\in (125T,128T)$ is given at the bottom: the vertical lines correspond to time instants at which the free surface is depicted in figure 25. As visible, the left side (orange solid) and the right side (black dash-point) signals are different and the intensity of the elevation switches between left and right after intervals of approximately $50\unicode{x2013}70$ oscillation periods. In figure 24(b), for example, the right side elevation is higher than the left side.

Figure 24. (a) Time histories of the surface elevations from the left probe $\eta _1$ and the right probe $\eta _{99}$ for the chaotic case $A=0.03 L$ and $T/T_1=0.5$. (b) Enlarged view of the time history for the probe $\eta _{99}$ within the time range $t\in [125T,128T]$. Left probe signal is drawn with the orange solid line, right probe with the black dash-point line. Vertical dashed lines correspond to time instants selected for the free-surface configurations of figure 25 and letters refer to related panels.

Figure 25. Chaotic case $A=0.03 L$ and $T/T_1=0.5$: free-surface configuration for 9 different time instants in two periods of the tank oscillations $t\in [125.75T, 127.75T]$. To improve the visibility $y^*$ is equal to $y$ for the red points, and equal to $y^*=y+0.15$ and $y^*=y+0.3$ for the green and blue points, respectively. Red points correspond to the configuration where $\eta _{99}$ is close to its maximum values of figure 24(b). Red arrows indicate the direction of the tank motion. The video of the simulation is available at Link Video N2.

The free surface shows different patterns when the tank is moving leftwards and rightwards. With respect to the interval highlighted in figure 24(a), when the free surface impinges on the left side of the tank and moves rightwards, a breaking wave is formed (see blue lines in figure 25). The same pattern is not formed when the free surface moves leftwards, as indicated by red lines in the same figure.

8. Sloshing dissipation

In the present section the slosh dissipation is addressed. As underlined in Colagrossi, Bouscasse & Marrone (Reference Colagrossi, Bouscasse and Marrone2015), the dissipation is linked to two main phenomena: (i) the wall boundary layer (WBL) near the tank walls and (ii) complex motions of the free surface leading to its fragmentation during breaking waves and wave impact events. These two mechanisms act with different intensity, depending on the attained regimes. Indeed, for the periodic monochromatic regime the motion of the free surface is rather mild and singular phenomena do not occur: as a consequence, the main dissipation comes from the WBL. Conversely, during the quasi-periodic regime, the large motion of the free surface, its fragmentation and the occurrence of water impact events make mechanism (ii) as the dominant one.

It is worth noting that the 2-D assumption does not take into account the viscous friction on the front and rear tank sides, hence, the WBL effects are reduced. In order to partially recover this effect, viscous damping corrections are often used in the literature, following the pioneering work of Keulegan (Reference Keulegan1959). These corrections are generally useful for sloshing flows in shallow water, as discussed in Bouscasse et al. (Reference Bouscasse, Antuono, Colagrossi and Lugni2013) and in Antuono et al. (Reference Antuono, Bouscasse, Colagrossi and Lugni2012a). Unfortunately, this is not the case for the present work, where the problem is in intermediate depth conditions and, for this reason, we preferred to avoid any correction to the fluid viscosity.

As introduced in § 3.3, the energy balance is given by three general terms

(8.1)\begin{equation} \mathcal{E}_M-\mathcal{W}_{NF}=\mathcal{E}_{diss}, \end{equation}

where $\mathcal {E}_M$ is the mechanical energy of the fluid in the Ni-FoR moving with the tank, $\mathcal {W}_{NF}$ is the work performed by the non-inertial forces and $\mathcal {E}_{diss}$ is the energy dissipated by the fluid.

Figure 26(a) shows the time behaviour of the three energy terms for the test case $A=0.03L$ and $T=0.867T_1$ (tripling-period scenario) evaluated by the $\delta$-LES-SPH model at resolution $N=H/\Delta r=200$. The energies in (8.1) are made non-dimensional with the reference kinetic energy

(8.2)\begin{equation} \Delta \mathcal{E}= \displaystyle{\frac{1}{2}} m_l \left( 2{\rm \pi} \displaystyle{\frac{A}{T}} \right)^2 \displaystyle{\frac{t_{end}}{T}}, \end{equation}

that is the maximum kinetic energy of the liquid (imagining it frozen inside the tank) multiplied by the total number of periods. From this figure it is clear that the mechanical energy needs approximately $N_{P0}=80$ oscillation periods before reaching a constant condition. Similarly, the external work $\mathcal {W}_{NF}$ and the dissipated energy $\mathcal {E}_{diss}$ time trends become almost linear. Within this regime the almost constant time rate of $\mathcal {E}_{diss}$ is defined as the dissipated power $\bar {\mathcal {P}}_{diss}$, which is directly related to the specific case investigated (see (3.14) in § 8).

Figure 26. Test case $A=0.03 L$ and $T/T_1=0.867$. (a) Time evolution of mechanical energy $\mathcal {E}_M$ (in Ni-FoR), work performed by the non-inertial forces $\mathcal {W}_{NF}$ and dissipated energy $\mathcal {E}_{diss}$. Results refer to $\delta$-LES-SPH simulations at finest spatial resolution $N=H/\Delta r=200$. (b) Time evolution of the energy dissipated by the fluid $\mathcal {E}_{diss}$ at three spatial resolutions $N=50$, 100 and 200. The energies are made non-dimensional with the reference energy in (8.2).

Panel (b) of the same figure depicts $\mathcal {E}_{diss}$ at three different spatial resolutions. The time trend of the dissipated energy clarifies that, at the finer resolution, the curve slope is convergent; moreover, it is also clear that $\mathcal {E}_{diss}$ is roughly well captured even at the lowest resolution, thus indicating that the present numerical approach is very reliable when energy evaluations are carried out.

By extending the discussion to the whole frequency spectrum and to both sloshing amplitudes, figure 27 shows $\bar {\mathcal {P}}_{diss}$ as a function of the oscillation periods and for both the motion amplitudes $A=0.01L$ and $0.03L$. The dissipated power is made non-dimensional with the reference power $\Delta \mathcal {P}=\Delta \mathcal {E}/t_{end}$.

Figure 27. Time rate of the dissipated energy distribution with frequency for $A=0.01 L$ and $A=0.03 L$. Results refer to $\delta$-LES-SPH simulations at the finest spatial resolution $N=H/\Delta r=200$.

The two curves, in logarithmic scale on the vertical axis, follow a similar trend as the WEFDs, discussed in §§ 4.1 and 4.2. It is worth noting that, where the tripling-period and the doubling-frequency modes occur, clearly visible peaks are found. The higher dissipation is linked to the presence of breaking waves and of wave impacts against the vertical walls, developing during the tripling-period and doubling-frequency regimes, respectively. The breaking wave phenomena are also responsible for the increase (in absolute value) of the dissipated power at the lowest excitation periods, where, at $A=0.03L$, chaotic regimes occur.

Finally, the time histories of the velocity of the fluid centre of mass $\dot {x}_G(t)$ and of the tank acceleration $a_{tank}(t)$ are depicted in figure 28 for $A=0.03L$ at $T=0.791 T_1$ and $T=1.098 T_1$. The former is one of the cases where the slosh dissipation is very low, hence, $\dot {x}_G(t)$ and $a_{tank}(t)$ are close to quadrature. Conversely, at the excitation period $T=1.098 T_1$ the maximum dissipation is attained and $\dot {x}_G(t)$ and $a_{tank}(t)$ are in anti-phase. Those phase lag conditions are in agreement with the formula (3.14), which links the average slosh dissipation rate $\bar {\mathcal {P}}_{diss}$ with $\dot {x}_G(t)$ and $a_{tank}(t)$.

Figure 28. Test cases $A=0.03L$, $T=0.791 T_1$ (a), $T=1.098 T_1$ (b) showing time histories of the fluid centre of mass velocity $\dot {x}_G(t)$ and the tank acceleration $a_{tank}(t)$.