5.1. The CPU costs

All simulations were performed on GLiCID, the Centrale Nantes supercomputing facility equipped with 96-core AMD EPYC Genoa processors (9424F series). To assess numerical convergence, spatial resolutions of $N = 22, 33, 50$ and $75$ were adopted for the three-dimensional simulations. At the highest resolution ( $N = 75$ , corresponding to $N_{\textit{part}} = 2.8 \times 10^{6}$ particles), the run required approximately 4 days of wall-clock time on 288 cores, with $N_{\textit{iter}} = 3.5 \times 10^{5}$ time steps.

Two-dimensional simulations were executed on a single GLiCID node (96 cores). For the two-dimensional convergence analysis, five resolutions were considered: $N = 50, 100, 200, 400$ and $800$ . At the finest resolution ( $N = 800$ , $N_{\textit{part}} = 2.14 \times 10^{6}$ ), the simulation required approximately 20 days of wall-clock time and $N_{\textit{iter}} = 1.65 \times 10^{6}$ iterations.

A three-dimensional simulation at $N = 800$ would correspond to approximately $3.41 \times 10^{9}$ particles and nearly four years of computing time on the full GLiCID facility (48 nodes, 4608 cores), which remains computationally prohibitive with current resources. Future work will explore adaptive-resolution strategies and GPU acceleration to mitigate these constraints.

Due to the chaotic nature of the flow, at least ten independent simulations were performed at each resolution to obtain reliable ensemble averages of the sloshing dissipation $\mathcal{Q}_{\textit{diss}}$ . As discussed in Marrone et al. (Reference Marrone, Colagrossi, Calderon-Sanchez and Martinez-Carrascal2021a ), the uncertainty decreases with increasing particle number; for $N = 50$ , the uncertainty is approximately $1\,\%$ . Therefore, only a single simulation was performed at the highest resolutions for the two-dimensional (2-D) ( $N = 800$ ), RW-3-D ( $N = 100$ ) and 3-D ( $N = 75$ ) cases.

In total, 31 simulations were carried out in three dimensions, 11 in RW-3-D and 82 in two dimensions (including cases with and without surface tension), for a total of 124 runs. These figures highlight the severe computational cost of high-resolution simulations, which currently limits the feasibility of fully resolved 3-D runs at extreme resolutions.

5.2. Three-dimensional results: global loads and energy dissipation

The 3-D simulations provide insight into the global hydrodynamic loads acting on the tank and the associated energy-dissipation mechanisms. Owing to computational constraints, the maximum spatial resolution was limited to $N = 75$ .

Figure 4.

Free surface at the first roof impact in the 3-D simulation. (a) three-dimensional view; (b) lateral and top views.

Figure 5.

Comparison between experimental data and the 3-D $\delta$ -LES–SPH simulation of the vertical force (top) and its dynamic component $F^{\textit{dyn}}_{\!y}$ (see (5.2)) acting on the tank. For the simulation, the spatial resolution $N = 75$ was adopted.

The first roof impact, occurring at $t = 1.08T$ , is shown in figure 4, which highlights the fully 3-D character of the free-surface dynamics. During the initial upward acceleration of the tank, a Rayleigh–Taylor instability develops and is modulated by surface waves generated by the collapse of the menisci along the lateral walls. These waves propagate inward from the sides towards the centre of the domain. As the flow evolves, the interface becomes increasingly complex, exhibiting a highly 3-D topology and intricate interactions among vortex tubes, which lead to a dynamically rich flow field.

The comparison between the 3-D $\delta$ -LES–SPH results and the experimental measurements of the vertical force $F_{\!y}$ acting on the tank is shown in figure 5(a). As previously observed in Marrone et al. (Reference Marrone, Colagrossi, Calderon-Sanchez and Martinez-Carrascal2021a ,Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez b ), the dominant contribution to $F_{\!y}$ is due to the liquid inertia, $M_L\,(g + a_{\textit{tank}})$ . Therefore, a more meaningful comparison between numerical and experimental data is obtained by introducing the dynamic force component

(5.2) \begin{equation} F^{\textit{dyn}}_{\!y}(t) = F_{\!y}(t) + M_L\,\left [g + a_{\textit{tank}}(t)\right ]\!. \end{equation}

As discussed in Malan et al. (Reference Malan, Pilloton, Colagrossi and Malan2022), this quantity corresponds to the inertial response of the liquid in the NiFoR, i.e. $F^{\textit{dyn}}_{\!y} = -M_L\,\ddot {y}_G$ , and is therefore non-zero only when sloshing develops inside the tank. The comparison between numerical and experimental time histories of $F^{\textit{dyn}}_{\!y}$ (figure 5 b) reveals discrepancies that are more evident than for the total force. Nevertheless, the overall agreement remains satisfactory, especially considering that the simulations were performed using a single-phase formulation and without explicit modelling of air cushioning.

5.2.1. Vortex identification and energy dissipation

Figure 6.

A $\delta$ -LES–SPH 3-D simulation. Rendering of the free-surface configuration (a) and iso-surface of the $Q$ -criterion, $Q=50\,\varOmega ^2$ (b), during the first bottom liquid impact ( $t=1.65T$ ).

Figure 6. Long description

Panel A: A 3-D simulation rendering of the free-surface configuration during the first bottom liquid impact. The simulation shows a rectangular container filled with a turquoise fluid, with the fluid surface exhibiting complex, turbulent behavior. The fluid appears to be in motion, creating waves and splashes within the container. Panel B: A 3-D simulation rendering of the iso-surface of the Q-criterion during the first bottom liquid impact. The simulation shows a rectangular container filled with a red fluid, with the iso-surface highlighting the vortical structures within the fluid. The fluid appears to be in turbulent motion, with intricate patterns and swirls visible within the container.

The well-known $Q$ -criterion (Hunt et al. Reference Hunt, Wray and Moin1988)

(5.3) \begin{equation} Q := \frac {1}{4}\left (\|\boldsymbol{\nabla }\times \boldsymbol{u}\|^2 - \|\mathbb{D}\|^2\right )\!, \end{equation}

is adopted to identify vortical structures within the flow. Here, $\|\mathbb{D}\|$ is the rescaled Frobenius norm of the strain-rate tensor already introduced in § 2.4. Positive values of $Q$ indicate regions where rotation dominates over strain, thus revealing coherent vortical structures.

The first water impact against the tank bottom occurs at $t = 1.65T$ (figure 6 a). The right panel of the same figure shows iso- $Q$ surfaces illustrating the complex vortex tubes generated during the sloshing process. The threshold $Q = 50\,\varOmega ^2$ was selected to highlight the dominant vortex structures while filtering out unresolved small scales. This specific time instant was chosen because, during water-impact events, the liquid is more confined and the free surface is less fragmented, making vortices more clearly visible. It is worth noting that, at the adopted resolution ( $N = 75$ ), the extracted iso- $Q$ surfaces remain sensitive to the chosen threshold, preventing an objective definition of turbulent eddy size.

The non-dimensional dissipated energy $\mathcal{Q}_{\textit{diss}}/\varDelta \mathcal{E}$ , obtained by increasing the spatial resolution from $N=22$ to $N=75$ , is reported in figure 7. Owing to the chaotic nature of the flow and following the analysis suggested in Marrone et al. (Reference Marrone, Colagrossi, Calderon-Sanchez and Martinez-Carrascal2021a , Reference Marrone, Colagrossi, Gambioli and González-Gutiérrezb ), ten independent simulations were performed at each resolution, except for $N=75$ , for which only one realisation was carried out because of the computational cost. The curves indicate progressive convergence towards the experimental reference, with a convergence rate of approximately 1.6. Despite the relatively coarse resolution and the single-phase model, the numerical predictions agree well with the experiments, confirming the robustness of the proposed framework for global quantities.

Although global integral quantities such as energy dissipation converge satisfactorily at moderate resolutions, the detailed characterisation of local flow features – including vortex topology and surface-tension effects – requires significantly higher resolution. This limitation is evident in figure 8, where the 3-D simulation at $N=75$ is compared with a 2-D simulation at $N=800$ , clearly showing the enhanced small-scale resolution that can be achieved in two dimensions.

Figure 7.

Convergence analysis of the dissipated energy $\mathcal{Q}_{\textit{diss}}$ in 3-D simulations performed with the $\delta$ -LES–SPH model. (a) time histories of slosh dissipation for four spatial resolutions and experimental data from Martinez-Carrascal & Gonzalez (Reference Martinez-Carrascal and Gonzalez2022). (b) relative error versus resolution $N$ .

Figure 8.

A $\delta$ -LES–SPH contour plot of the $Q$ -criterion at $t=1.65T$ . (a) cross-section at $z=W/6$ of the 3-D simulation at $N=75$ . (b) Two-dimensional simulation at $N=800$ with surface tension.

Figure 8. Long description

Panel A: A heat map representing the Q-criterion distribution in a cross-section of a 3-D simulation at a specific time. The x-axis ranges from 0 to 0.1 meters, and the y-axis ranges from 0 to 0.06 meters. The color scale ranges from 0 to 50, with blue indicating lower values and red indicating higher values. The heat map shows a varied distribution of the Q-criterion, with some regions having higher intensity (red) and others having lower intensity (blue). Panel B: A heat map representing the Q-criterion distribution in a two-dimensional simulation at a specific time with surface tension. The x-axis ranges from 0 to 0.1 meters, and the y-axis ranges from 0 to 0.06 meters. The color scale ranges from 0 to 50, with blue indicating lower values and red indicating higher values. The heat map shows a more dispersed distribution of the Q-criterion compared to Panel A, with noticeable clusters of higher intensity (red) and lower intensity (blue) regions.

5.3. Reduced-width 3-D simulations: direct evidence on the role of the spanwise dimension

The 2-D analysis carried out in the previous sections rests on the physical hypothesis, supported by the literature on bounded violent free-surface flows (Lubin et al. Reference Lubin, Vincent, Abadie and Caltagirone2006; Marrone et al. Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez2021b ), that, in this class of problems, the dominant mechanism for the generation and destruction of coherent vortical structures is set by jet and droplet impacts against the walls and by fragmentation and reconnection of the interface, rather than by intrinsically 3-D processes associated with vortex stretching.

Figure 9.

An $\delta$ -LES–SPH comparison of the vortex-tube configurations for the reduced-width 3-D simulation (top) and the full 3-D simulation (bottom). Free-surface configuration and iso-surface of the $Q$ -criterion, $Q=100\,\varOmega ^2$ , during the second (a,c) and third (b,d) bottom liquid impacts ( $t=2.58T$ and $t=3.63T$ ). Both simulations are performed with the spatial resolution $N=100$ .

To provide direct evidence in support of this hypothesis, a dedicated set of RW-3-D simulations was performed by reducing the spanwise length by a factor of five, namely $W_z = W/5 = 12$ mm. Free-slip boundary conditions are imposed on the two spanwise walls, $z=0$ and $z=W_z$ , to avoid interactions between the boundary layers on those solid boundaries in the narrow domain. The adopted spatial resolutions are $N = 50$ and $N = 100$ (see table 1). In the RW-3-D configuration, the spanwise direction is fully resolved, and vortex stretching is therefore allowed wherever physically active. Although the reduction of the spanwise extent makes it possible to reach a resolution of $N=100$ , surface-tension effects still cannot be adequately resolved and are therefore not included in the RW-3-D framework.

Figure 9(a,b) shows the free-surface configurations together with an iso-surface of the $Q$ -criterion ( $Q = 100\,\varOmega ^2$ ) at two representative instants of the RW-3-D simulation at $N = 100$ , corresponding to the second and third bottom liquid impacts ( $t = 2.58\,T$ and $t = 3.63\,T$ ). At both instants, the flow downstream of the impact is densely populated by 3-D vortical structures that span the entire reduced width $W_z$ and exhibit clear signatures of spanwise stretching and bending. The reduction of the spanwise extent therefore does not suppress the three-dimensional character of the small-scale dynamics: vortex stretching and the associated production of 3-D vorticity are admitted by the discretisation and develop spontaneously wherever they are physically active.

To make the comparison of the vortical structures as direct as possible, the full-width 3-D simulation used in figure 9 was also performed at the same spatial resolution as the RW-3-D case, namely $N=100$ . Owing to the very high computational cost, this additional full-width 3-D run was advanced only up to $t=3.63T$ , which is sufficient to compare the vortex-tube configurations at the second and third bottom impacts.

Figure 9(c,d) depicts the corresponding free-surface configurations and iso-surfaces of the $Q$ -criterion ( $Q = 100\,\varOmega ^2$ ) at the same time instants for the full 3-D simulations. The vortex tubes appear longer and more intricate because of the larger spanwise dimension. Consequently, vortex stretching is more pronounced, while free-surface fragmentation is less energetic, owing to the larger space available for the free surface to deform.

Despite these clear differences in the vortex distributions obtained in the RW-3-D and full 3-D simulations, their effect on slosh dissipation appears negligible. This is shown in figure 10(a), which reports the time histories of the non-dimensional slosh dissipation $\mathcal{Q}_{\textit{diss}}/\varDelta \mathcal{E}$ for the RW-3-D simulations at $N = 50$ and $N = 100$ , together with the full 3-D simulation at $N = 75$ and the experimental reference. The four curves are mutually consistent throughout the entire simulation window: the two RW-3-D resolutions are already converged with respect to each other at $N = 100$ , and both agree closely with the full 3-D result and with the experimental data, with deviations of the order of a few per cent throughout the sloshing history. The reduction of the spanwise extent by a factor of five with respect to the original tank therefore leaves the dimensionless global energy dissipation unchanged within the experimental uncertainty. It is thus possible to conclude that, for this violent free-surface flow, despite the presence of a genuinely 3-D dynamics, the integral quantities of the flow are essentially unaffected by the spanwise reduction.

Figure 10.

A $\delta$ -LES–SPH reduced-width 3-D simulation. (a) time histories of slosh dissipation for the reduced-width 3-D simulations at resolutions $N=50$ and $N=100$ , compared with the full 3-D simulation ( $N=75$ ) and the experimental data. (b) contour plot of the $Q$ -criterion in the cross-section at $z=W/8$ of the reduced-width 3-D simulation at $N=100$ at the first bottom impact ( $t=1.65T$ ).

Figure 10. Long description

The image contains two panels labeled (a) and (b). Panel (a) is a line graph showing time histories of slosh dissipation. The x-axis represents normalized time (t/T) and the y-axis represents normalized dissipation (Qdiss/Δε). The graph includes four data series: RW_3-D N = 50 (green dashed line), RW_3-D N = 100 (blue dashed line), 3-D N = 75 (orange dashed line), and experimental data (black solid line). The data series show a general trend of decreasing dissipation over time, with the experimental data closely followed by the 3-D N = 75 simulation. Panel (b) is a contour plot showing the Q-criterion in a cross-section at z = W/8 of the reduced-width 3-D simulation at t = 1.65T. The x-axis represents the x-coordinate in meters, and the y-axis represents the y-coordinate in meters. The contour plot uses a color scale to represent the Q/Ω^2 values, ranging from blue (0) to red (50.0). The plot shows a complex turbulent flow pattern with regions of high and low Q-criterion values.

Figure 10(b) shows a slice of the $Q$ -criterion at $z = W/8$ for the RW-3-D simulation at $N = 100$ , taken at the first bottom impact ( $t = 1.65\,T$ ). Comparison with the fields shown in figure 8 reveals that the vortical structures span a wider range of scales. In particular, some of the largest structures visible in this slice may result from a geometrical sampling effect: in the reduced-width configuration, a transverse cut has a higher probability of intersecting longitudinal vortices developing along the $x$ -axis than in the full-width 3-D domain.

The observations discussed in this section have a direct physical interpretation that supports the central thesis of this work. Although 3-D vortical structures and vortex stretching are clearly present in the RW-3-D simulation (figure 9), their presence does not translate into a measurable modification of the energy budget (figure 10). The bulk turbulent dynamics in the spanwise direction is therefore essentially a passive consequence of the in-plane dynamics, which is in turn set by the bounded geometry and by the breaking free surface. The 2-D simulations should thus not be regarded merely as an approximation of the 3-D case: in the present class of bounded violent free-surface flows, the two descriptions converge on the quantities that govern the engineering behaviour of the system. This is a property of the physical problem, not a numerical coincidence, and it has a direct positive consequence: the much lower computational cost of the 2-D framework grants access to in-plane resolutions ( $\Delta r = H/800$ at $N = 800$ ) that are currently out of reach in fully 3-D simulations of this kind of violent free-surface flow.

5.4. Two-dimensional results: turbulence analysis and vortex detection (case with $\textit{We}=\infty$ )

Having established in the previous subsections that the 3-D and RW-3-D simulations yield consistent global loads and energy dissipation, the analysis now turns to the 2-D framework.

As discussed above, in the present class of bounded violent free-surface flows, the 2-D description faithfully reproduces the integral quantities of interest, while granting access to a substantially finer in-plane resolution than is affordable in three dimensions. This higher resolution is exploited here to characterise in detail the turbulent scales, the genesis and lifespan of vortical structures and the spatial distribution of the eddies.

The 2-D analysis initially neglects surface-tension effects ( $\textit{We} = \infty$ ). However, as noted in Marrone et al. (Reference Marrone, Colagrossi, Calderon-Sanchez and Martinez-Carrascal2021a ,Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez b ), the presence of menisci at the intersection between the lateral walls and the liquid surface plays a crucial role in triggering the initial disturbance of the otherwise flat free surface. After the first upward motion of the tank, during the first downward motion, i.e. for $t\gt T$ (see figure 2), a Rayleigh–Taylor instability develops, initiated by the meniscus geometry. To account for these effects, the menisci are numerically replicated in the initial condition by adding a small fluid triangle near the vertical boundaries, forming a contact angle consistent with water surface tension.

Figure 11.

A 2-D $\delta$ -LES–SPH vorticity field at six different time instants: (a) tank near maximum vertical position, (b) first roof impact, (c) first bottom impact, (d) third roof impact, (e) third bottom impact and (f) twelfth bottom impact. Link to Video 1.

Figure 12.

Comparison between experimental data and the 2-D $\delta$ -LES–SPH simulation of the vertical force acting on the tank ( $N = 800$ ).

Figure 13.

Comparison between experimental data and the $\delta$ -LES–SPH simulation of the vertical dynamic component $F^{\textit{dyn}}_{\!y}$ (see (5.2)) for the 2-D case at $N = 800$ .

Figure 11 shows six representative instants of the 2-D simulation at the finest spatial resolution ( $N = 800$ ), reported in the tank reference frame. Particles are coloured according to their vorticity. As in the 3-D case, the flow first exhibits a Rayleigh–Taylor instability during the initial downward motion, modulated by surface waves generated by the collapse of the menisci along the lateral walls. These waves propagate towards the centre of the domain, promoting progressive fragmentation of the free surface and inducing repeated liquid impacts against the horizontal walls of the tank. A video of the simulation is available through the link provided in the caption of figure 11.

Figures 12 and 13 report the time histories of the vertical force $F_{\!y}$ and its dynamic component $F_{\!y}^{\textit{dyn}}$ , compared with experimental measurements. The vertical dashed lines indicate the six instants corresponding to the snapshots in figure 11. The agreement between numerical and experimental data is satisfactory, confirming that the $\delta$ -LES–SPH model captures the global sloshing dynamics even in the absence of surface tension.

The 2-D framework allows for much finer spatial resolution than the 3-D case, enabling a more detailed characterisation of turbulence and the vortex dynamics. In the next section, we analyse the spectral content of the flow and quantify the distribution of kinetic energy across scales.

5.4.1. Analysis of the turbulence scales

The analysis of turbulence scales provides insight into the resolution requirements for accurately capturing the dynamics of violent-sloshing flows. This is particularly relevant for assessing the adequacy of the $\delta$ -LES–SPH approach compared with direct numerical simulation (DNS).

Figure 14.

Two-dimensional $\delta$ -LES–SPH. Top: time histories of kinetic and potential energies in the NiFoR. Bottom: time histories of vertical tank displacement and energy dissipation rate $\mathcal{P}_{\textit{diss}}$ .

Figure 15.

Comparison of kinetic energy and dissipation rate between 2-D ( $N=800$ ) and 3-D ( $N=75$ ) simulations.

Figure 15. Long description

Two line graphs compare kinetic energy and dissipation rate between 2-D and 3-D simulations. Panel A: The line graph shows the kinetic energy over time. The x-axis represents time normalized by a period (t/T) ranging from 0 to 10. The y-axis represents the normalized kinetic energy (ε_K/Δε) ranging from 0 to 1.2. The black solid line represents 2-D simulations with 800 particles, while the orange dotted line represents 3-D simulations with 75 particles. The graph shows fluctuating patterns for both simulations, with the 3-D simulation exhibiting higher peaks and more frequent fluctuations compared to the 2-D simulation. Panel B: The line graph shows the dissipation rate over time. The x-axis represents time normalized by a period (t/T) ranging from 0 to 10. The y-axis represents the normalized dissipation rate (-P_diss/ΔP) ranging from 0 to 10. The green solid line represents 2-D simulations with 800 particles, while the orange dotted line represents 3-D simulations with 75 particles. The graph shows frequent spikes in dissipation rate for both simulations, with the 3-D simulation displaying more pronounced and frequent spikes compared to the 2-D simulation.

The analysis is performed at $t = 3.52T$ (frame ‘e’ in figure 11). Figure 14 reports the dimensionless turbulent kinetic energy $\mathcal{E}$ and the rate of energy dissipation $\mathcal{P}_{\textit{diss}}$ . The turbulent kinetic energy and dissipation per unit mass, $k$ and $\epsilon$ , are non-dimensionalised as

(5.4) \begin{align} \tilde {k}=\frac {k}{gH}= \left (\frac {\mathcal{E}_K}{\varDelta \mathcal{E}}\right )\left (\frac {2A}{H}\right )\!, \qquad \tilde {\epsilon }=\frac {\epsilon }{g\sqrt {gH}}= \left (\frac {-\mathcal{P}_{\textit{diss}}}{\varDelta \mathcal{P}}\right ) \left (\frac {2A\,f_0}{\sqrt {gH}}\right )\!. \end{align}

From figure 14, at $t = 3.52T$ , it follows that $\mathcal{E}_K / \varDelta \mathcal{E} \simeq 0.72$ and $\mathcal{P}_{\textit{diss}} / \varDelta \mathcal{P} \simeq -4.1$ . Consequently, $\tilde {k} \simeq 2.7$ and $\tilde {\epsilon } \simeq 5.6$ .

5.4.1.1. Integral scale

Following Marrone et al. (Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez2021b ), the scale of the large eddies is

(5.5) \begin{equation} \ell =\frac {k^{3/2}}{\epsilon }=\frac {\tilde {k}^{3/2}}{\tilde {\epsilon }}\,H\simeq 0.8H, \end{equation}

indicating that the largest vortices have a size comparable to the filling height $H$ .

5.4.1.2. Taylor microscale

The Taylor microscale $\lambda$ marks the boundary of the inertial subrange and is given by

(5.6) \begin{equation} \lambda =\sqrt {10\frac {\nu k}{\epsilon }}= H\sqrt {10\frac {\tilde {\nu }\tilde {k}}{\tilde {\epsilon }}}\simeq 0.017H=\frac {H}{57}, \end{equation}

where $\tilde {\nu }=\nu /(H\sqrt {gH})$ . The adopted resolutions ( $N=800$ in two dimensions and $N=75$ in three dimensions) are both larger than $H/\lambda =57$ , ensuring that the primary turbulent eddies are resolved.

Figure 16.

Two-dimensional $\delta$ -LES–SPH simulation at $t = 3.52T$ , resolution $N = 800$ . (a) vorticity field. (b) eddy-viscosity field. Note that $\mu _T$ is approximately half the value reported in Michel et al. (Reference Michel, Durante, Colagrossi and Marrone2022) because of the doubled resolution adopted in this work.

5.4.1.3. Kolmogorov microscale

The smallest scale is

(5.7) \begin{equation} \eta =\left (\frac {\nu ^3}{\epsilon }\right )^{1/4}= H\left (\frac {\tilde {\nu }^3}{\tilde {\epsilon }}\right )^{1/4}\simeq \frac {H}{2200}. \end{equation}

Resolving $\eta$ would require DNS with $N\approx 2200$ , which is computationally prohibitive in three dimensions. Only the 2-D simulation at $N=800$ approaches this limit.

The comparison in figure 15 shows that the time histories of $\mathcal{E}$ and $\mathcal{P}_{\textit{diss}}$ are generally consistent between two and three dimensions, with slightly higher kinetic energy and delayed peaks in the 3-D case. These results confirm that the adopted resolutions are sufficient to capture the large-scale dynamics, while the smallest scales remain unresolved in three dimensions.

Figure 16 shows the vorticity field (a) and the corresponding eddy viscosity (b) at $t = 3.52T$ for the 2-D simulation. The ratio between turbulent and molecular viscosity is of order one, confirming that the inertial scales are well resolved. These findings highlight the capability of the $\delta$ -LES–SPH approach to reproduce the main turbulent scales in the 2-D framework, while this remains impractical for 3-D violent-sloshing flows. The next section analyses the 2-D spatial distribution of the different eddies to further characterise the turbulence regime under investigation.

5.4.2. Vortex detection procedure

Figure 17 displays the contour plot of the $Q$ -criterion for the 2-D simulation at the analysed time instant $t = 3.52T$ (see (5.3)).

Figure 17.

(a) Two-dimensional $\delta$ -LES–SPH contour plot of the $Q$ -criterion at time $t = 3.52T$ . (b) enlarged view.

Figure 18.

(a) points of the fluid domain where the $Q$ -criterion exceeds $50\varOmega ^2$ . (b) vortex-size distribution normalised by the kernel support radius $l = 4\Delta r$ .

It is worth noting that the value of $Q$ inside the vortices is of order $\mathcal{O}(10^4)\,\varOmega ^2$ . Consequently, the contour plot shown in figure 17 remains essentially unchanged when the threshold value of $Q$ is varied within the range $0.1$ – $100\,\varOmega ^2$ . This invariance is a significant result, as it enables an objective definition of eddy size.

However, because of the lower spatial resolution, this behaviour is not observed in the 3-D framework, where the topology of the iso- $Q$ surfaces is more sensitive to the chosen $Q$ threshold.

The observed difference between 2-D and 3-D vortex detection is primarily attributed to the spatial resolution of the numerical simulations. In the 2-D case, the vortices are sharply resolved and the $Q$ -field exhibits steep gradients around the vortex cores. As a result, the topology of the iso- $Q$ contours remains nearly invariant when the visualisation threshold is varied over several orders of magnitude. Conversely, in the 3-D framework, the lower spatial resolution leads to a more diffuse $Q$ distribution, where rotational and strain-dominated regions are less clearly separated. This causes the iso- $Q$ surfaces to be highly sensitive to the chosen threshold, significantly affecting the apparent size and shape of the detected vortices. Such sensitivity highlights the importance of resolution in accurately capturing coherent structures in violent-sloshing flows.

Figure 18(b) shows the distribution of the vortical size $d_Q$ , normalised by the kernel support radius $l = 4\Delta r$ . A total of 2000 eddies are detected. The highest bar corresponds to small eddies whose size is smaller than $l$ ; these cannot be resolved with the current spatial resolution ( $N = 800$ ), and their effects on vorticity dissipation are accounted for by the subgrid model. Conversely, the remaining larger vortices are fully resolved.

5.4.3. Genesis, merging and break-up of the vortex structures

In the absence of a free surface, an initial array of $64$ Taylor–Green vortices in a doubly periodic 2-D domain (figure 19) progressively evolves through a sequence of pairwise merging events until only two large coherent structures survive. The dimensionless time in figure 19 is defined as $t\,U/L$ , where $U$ is the maximum amplitude of the velocity field in the initial condition and $L$ is the size of the square bi-periodic domain.

Figure 19.

Time evolution of the vorticity field of 2-D Taylor–Green vortices (adapted from Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021a )). The initial condition (top-left panel) consists of $64$ vortices arranged on an $8 \times 8$ array in a doubly periodic domain at $\textit{Re}=10^6$ . The six panels show snapshots at dimensionless times $tU/L$ equal to 0.0, 1.35 and 2.05 ((a–c) from left to right), and 3.30, 5.80 and 11.3 ((d–f), from left to right). Through a sequence of pairwise merging events, the system evolves towards two coherent vortical structures, as predicted by the theory of the 2-D inverse cascade. Link to Video 2.

Figure 20.

Genesis of vortex structures from jet and droplet impacts and their subsequent merging. The small yellow circle indicates the Taylor-microscale size.

Figure 20. Long description

The image contains nine panels, each depicting a line graph with contour plots. The x-axis represents the horizontal distance in meters (m), ranging from 0 to 0.020 m. The y-axis represents the vertical distance in meters (m), ranging from 0 to 0.035 m. Each panel shows the evolution of vortex structures at different time points, indicated by the variable t. The small yellow circle in each panel indicates the Taylor-microscale size. The contour plots use a color scheme to represent the spatial gradients, with red and blue indicating different gradient directions or magnitudes. Panel A shows the initial state at t = 2.89T. Panel B shows the state at t = 2.97T, where vortex structures begin to form. Panel C shows the state at t = 3.01T, with more defined vortex structures. Panel D shows the state at t = 3.05T, where vortex structures become more complex. Panel E shows the state at t = 3.07T, with further development of vortex structures. Panel F shows the state at t = 3.09T, where vortex structures start merging. Panel G shows the state at t = 3.10T, with continued merging of vortex structures. Panel H shows the state at t = 3.11T, where merging is more pronounced. Panel I shows the state at t = 3.12T, with advanced merging and complex vortex patterns.

The spontaneous upscale transfer is a direct manifestation of the inverse cascade, a distinctive feature of 2-D turbulence in the absence of boundaries. In 3-D flows, vortex stretching redirects energy towards smaller scales and prevents this long-lived coalescence. This reference case is reproduced here because it represents the baseline behaviour that the present sloshing flow would exhibit if the free-surface dynamics did not actively interfere with vortex evolution. As shown in figures 20 and 21, the interaction with the fragmented interface repeatedly disrupts the coalescence process, so that the quasi-stationary ‘two-vortex’ state of figure 19 is never attained in the present violent-sloshing regime. This contrast motivates the interpretation of the free surface as an active ingredient of the cascade dynamics.

This inverse energy cascade is also observed in the present 2-D sloshing simulations. Figure 20 shows the genesis and merging of vortex structures at nine time instants following the first roof impact. At $t = 2.89T$ , the vorticity field is primarily associated with wall boundary layers. At $t = 2.97T$ , water jets and droplets generated during the initial roof impact fall towards the underlying liquid. Their impact creates vortex structures that begin to convect within the fluid domain. This phenomenon has been extensively discussed in the SPH literature (see, e.g. Trivellato et al. Reference Trivellato, Bertolazzi and Colagrossi2006; Landrini et al. Reference Landrini, Colagrossi, Greco and Tulin2012; Antuono et al. Reference Antuono, Colagrossi, Le Touzé and Monaghan2013; Colagrossi et al. Reference Colagrossi, Bouscasse and Marrone2015).

Initially, the vortices in the top row of figure 20 have a size comparable to the Taylor microscale (marked with yellow circles). The process continues until $t = 3.05T$ , when a large cavity forms and begins to collapse. The collapse is completed at $t = 3.09T$ , producing an elongated vortex structure near the left tank wall. These vortices, much larger than the Taylor microscale, merge into even larger structures (see the lower panels of figure 20). This sequence explains how eddies are generated in this violent-sloshing flow, populating the eddy-size distributions discussed in the previous section.

Figure 21.

Tracking of a single vortex from the tank bottom to the free surface, where interface break-up disrupts the vortex and interrupts the inverse energy cascade.

Figure 21 illustrates the opposite process: large vortices breaking up upon reaching the free surface. This phenomenon has been analysed in detail by Chiron et al. (Reference Chiron, Marrone, Di Mascio and Le Touzé2018), who simulated the interaction of a dipole with a free surface and showed its fragmentation into smaller vortices. Similarly, Zhang et al. (Reference Zhang, Zhang, Zhang, Long, Han, Liu, Ohl and Li2025) revealed the mechanism driving the disruption of a vortex ring by a free surface. Although air plays a minor role in the present case, the underlying mechanism is consistent with the vortex break-up observed here.

An important observation is that the break-up of vortices upon reaching the free surface effectively limits the inverse energy cascade that characterises planar 2-D turbulence. In classical 2-D flows, energy is transferred from smaller to larger scales, leading to the formation of increasingly coherent structures. However, in violent-sloshing conditions, the interaction with the free surface introduces a strong dissipative mechanism: large vortices are disrupted and fragmented into smaller structures, interrupting the upscale energy transfer.

This process reduces the persistence of large-scale vortices and constrains the energy-transfer pathways that would otherwise distinguish 2-D and 3-D turbulence. Consequently, the free-surface interaction acts as a key factor in reducing the practical gap between idealised 2-D turbulence and a realistic sloshing dynamics.

This observation acquires additional significance in light of the independent recent study by Eswara & Sandilya (Reference Eswara and Sandilya2026), who showed that the two cascade mechanisms suppressed in 2-D flow – vortex stretching and strain-rate self-amplification – vanish not merely quantitatively but identically in two dimensions, concluding that 2-D simulations of sloshing flow underpredict damping by two orders of magnitude relative to 3-D results. The present violent-sloshing regime offers a complementary perspective: although these two cascade mechanisms are indeed absent in 2-D formulations, the free-surface-induced vortex destruction documented here constrains the energy-transfer pathways in both frameworks, reducing the practical gap in total dissipated energy between 2-D and 3-D simulations under impact-dominated conditions. The two results are therefore not in contradiction but address different flow regimes: Eswara & Sandilya (Reference Eswara and Sandilya2026) focus on a gentle sloshing flow at moderate Reynolds numbers where the direct cascade is the dominant dissipation pathway, whereas the present work considers violent vertical excitation, where impact and fragmentation events dominate and the free surface itself regulates the lifespan of turbulent eddies.

5.4.4. Surface-tension effects on the vorticity field

In this section, the role of surface tension in the 2-D framework is analysed. Figure 22 depicts the time evolution of the sloshing flow at Weber number $\textit{We} = 2240$ , corresponding to the experimental conditions (see § 4).

The presence of surface tension significantly reduces the fragmentation of the free surface. At the first instant, $t = 1.01T$ , the jets rise less than in the case $\textit{We} = \infty$ , resulting in a smoother interface compared with the corresponding case without surface tension. As the flow evolves, the jets eventually reach the top, but the free surface remains considerably less jagged and fragmented.

Figure 22.

Two-dimensional $\delta$ -LES–SPH vorticity field at six different time instants for $\textit{We} = 2240$ : (a) tank near maximum vertical position, (b) first roof impact, (c) first bottom impact, (d) third roof impact, (e) third bottom impact and (f) twelfth bottom impact. Link to Video 3.

Figure 22. Long description

The image contains six panels of vorticity fields at different time instants for a tank impacting surfaces. Each panel is a contour plot with color gradients representing vorticity values. The x-axis and y-axis are labeled in meters. The color bar indicates vorticity values ranging from -50 to 50. Panel A shows the tank near its maximum vertical position. Panel B captures the first roof impact. Panel C depicts the first bottom impact. Panel D illustrates the third roof impact. Panel E shows the third bottom impact. Panel F represents the twelfth bottom impact. Each panel captures the vorticity field at a specific time instant, with the time indicated at the top of each panel. The color gradients and contour lines show the distribution and intensity of vorticity in the fluid flow around the tank at each impact stage.

To provide a quantitative assessment, the analysis of turbulent scales performed in § 5.4.1 is repeated here. The selected time instant is $t = 3.67T$ , with a delay of $0.15T$ compared with the case $\textit{We} = \infty$ , since surface tension delays water-impact events. At $t = 3.67T$ , corresponding to the third liquid impact against the tank bottom, suitable conditions for turbulence analysis are achieved. The dimensionless turbulent kinetic energy and dissipation rate are $\mathcal{E}_K/\varDelta \mathcal{E}= 0.53$ and $\mathcal{P}_{\textit{diss}}/\varDelta \mathcal{P} = 2$ , respectively, corresponding to turbulent kinetic energy and dissipation per unit mass of $\tilde {k} \simeq 2$ and $\tilde {\epsilon } \simeq 2.73$ . The estimated Taylor microscale and Kolmogorov scale are

(5.8) \begin{align} \lambda = H/50, \qquad \eta = H/2000, \qquad (\textit{We}=2240). \end{align}

As expected, these scales are larger than in the case $\textit{We} = \infty$ , although the differences remain limited because of the large Weber number considered.

A comparison of the vorticity fields shown in figures 11 and 22 clearly indicates that surface tension reduces the number of eddies. A useful physical scale for the interplay between surface tension and inertial forcing is the capillary length $\ell _c$ (Eggers & Villermaux Reference Eggers and Villermaux2008)

(5.9) \begin{align} \ell _c = \sqrt {\frac {\sigma }{\rho _0 g}} = H\frac {\textit{Fr}}{\sqrt {\textit{We}}} \simeq \frac {H}{11} \approx 2.7\,\mathrm{mm}, \end{align}

for water at room temperature, where $g$ is the gravitational acceleration, which marks the scale below which surface tension dominates over gravity in shaping the free surface. In the present violent-sloshing problem, however, the relevant inertial force is set by the imposed tank acceleration $a_{\textit{max}} \approx 8.8\,g$ (cf. § 4) rather than by gravity alone, so that the capillary length is reduced to

(5.10) \begin{align} \ell _c = \sqrt {\frac {\sigma }{\rho _0 a_{\textit{max}}}} = H\frac {\textit{Fr}}{\sqrt {\textit{We}}}\frac {1}{\sqrt {a_{\textit{max}}/g}} \simeq \frac {H}{33} \approx 0.92\,\mathrm{mm}. \end{align}

This length is well resolved by the numerical discretisation at the finest resolution considered, where $\Delta r = H/800$ , i.e. $\ell _c/\Delta r \approx 24$ , and it is comparable to the typical size of the coherent droplets and jets visible in figure 22.

During the flow evolution, the effective capillary length ranges between the above two values, i.e. $\ell _c \in (0.92,2.7)$ mm. This provides a physical rationalisation of the reduced fragmentation observed with surface tension: the capillary scale sets a lower bound on the size of the free-surface structures that can be generated by impact-driven fragmentation, and consequently on the size of the smallest vortices generated at the interface. Indeed, $\ell _c \in (1.52,4.55)\,\lambda$ .

Figure 23 shows the vorticity field and the corresponding eddy viscosity. Consistently with the reduced turbulence intensity shown by the comparison of the vorticity fields, the eddy-viscosity field exhibits lower values when surface-tension effects are included.

Figure 23.

(a) Two-dimensional $\delta$ -LES–SPH vorticity at maximum resolution $N = 800$ , $t = 3.67T$ . (b) corresponding eddy viscosity. Weber number $\textit{We} = 2240$ .

Figure 24.

Two-dimensional $\delta$ -LES–SPH contour plot of the $Q$ -criterion at $t = 3.67T$ for $\textit{We} = 2240$ . (a) full view. (b) enlargement of the highlighted region.

Figure 24. Long description

Panel A: A heat map showing vorticity fields with different levels of surface tension. The x-axis ranges from 0 to 0.1 meters, and the y-axis ranges from 0 to 0.06 meters. The color scale ranges from blue to red, indicating values from 0 to 50. The map shows a complex distribution of vorticity with various eddies and regions of high and low vorticity. Panel B: An enlarged view of a specific region from Panel A. The x-axis ranges from 0.060 to 0.068 meters, and the y-axis ranges from 0.012 to 0.016 meters. The color scale is the same as in Panel A. The enlarged view highlights the Taylor microscale, which is marked as 0.6 millimeters. The map shows detailed structures of vorticity within this specific region.

Figure 25.

(a) points of the fluid domain where the $Q$ -criterion exceeds $50\varOmega ^2$ for the case with surface tension ( $\textit{We}=2240$ ), at $t = 3.67\,T$ . (b) vortex-size distribution $d_Q$ normalised by the kernel support radius $l = 4\Delta r$ , comparing the case with surface tension ( $\textit{We}=2240$ , green) and without surface tension ( $\textit{We}=\infty$ , red, reproduced from figure 18 b). The arrow marks the Taylor microscale $\lambda /l$ .

Figure 24 shows the contour plot of $Q$ for the analysed time instant $t = 3.67T$ . The enlarged view on the right highlights that the large eddies are well resolved.

Following the vortex detection procedure described in § 5.4.2, figure 25 shows the distribution of the vortical size $d_Q$ , normalised by the kernel support radius $l = 4\Delta r$ . A total of 2211 eddies are detected. The direct comparison of the two histograms in figure 25(b) shows that the two distributions are nearly indistinguishable for the smallest eddies (with size $d_Q \lesssim 4\,l$ , i.e. comparable to or below the Taylor microscale), whereas they differ markedly for larger structures. When surface tension is included ( $\textit{We}=2240$ ), the population of vortices with size $d_Q \gtrsim 7\,l$ is substantially higher than in the $\textit{We}=\infty$ case, by up to an order of magnitude in the largest bins. This is consistent with the reduced fragmentation of the free surface promoted by surface tension, which favours the genesis and persistence of larger coherent structures rather than altering the population of the smallest, SGS eddies.

In summary, surface tension acts as a stabilising mechanism, reducing free-surface fragmentation and turbulence intensity while promoting the persistence of larger coherent vortical structures.

5.5. Slosh dissipation analysis in the 2-D framework

As outlined in § 3, the total energy dissipation can be decomposed into two main contributions: the numerical dissipation $\mathcal{Q}_{N}$ and the viscous dissipation $\mathcal{Q}_{V} + \mathcal{Q}_{V}^{\textit{turb}}$ . The numerical dissipation originates from the density-diffusion term $\mathcal{D}_i^{\,\rho }$ (see § 2.5), which becomes active when strong pressure waves are generated after fluid–fluid or fluid–wall impacts. As emphasised by Antuono et al. (Reference Antuono, Marrone, Colagrossi and Bouscasse2015), this component theoretically converges to zero as the spatial resolution increases, i.e. as $\Delta r \to 0$ . In practice, however, as shown later in this section, even at the highest resolution ( $N = 800$ ), $\mathcal{Q}_N$ remains not only non-negligible but comparable to the physical dissipative contributions, $\mathcal{Q}_{V} + \mathcal{Q}_{V}^{\textit{turb}}$ .

Following the analysis of Malan et al. (Reference Malan, Pilloton, Colagrossi and Malan2022) and the study by Marrone et al. (Reference Marrone, Colagrossi, Di Mascio and Le Touzé2015), when an incompressible solver is used, $\mathcal{Q}_N$ is replaced by the energy jump associated with the pressure-Poisson equation. Indeed, liquid impacts represent singular events for an incompressible medium, leading to a lack of strict mechanical-energy conservation, as underlined, for example, by Szymczak (Reference Szymczak1994b ).

In the following, the analysis focuses on the viscous dissipation $\mathcal{Q}_{V} + \mathcal{Q}_{V}^{\textit{turb}}$ , which represents the dominant contribution to the total energy loss. To better characterise the dissipation induced by turbulent eddies, the discussion is shifted from $\mathcal{Q}_{V} + \mathcal{Q}_{V}^{\textit{turb}}$ to $\mathcal{Q}_{\omega } + \mathcal{E}_{\textit{Def}}$ , where $\mathcal{Q}_{\omega }$ denotes the cumulative energy dissipation associated with enstrophy, and $\mathcal{E}_{\textit{Def}}$ refers to the cumulative dissipation caused by viscous deformation of the fluid domain due to free-surface motion.

In general, considering the low viscosity of water, the velocity field of sloshing flows can be regarded as nearly irrotational and divergence free, implying that $\mathcal{E}_{\textit{Def}}$ should dominate over $\mathcal{Q}_{\omega }$ . However, this assumption does not hold in the present study: the violent flow conditions lead to extensive free-surface fragmentation and to the development of a turbulent regime, making $\mathcal{Q}_{\omega }$ comparable to $\mathcal{E}_{\textit{Def}}$ .

More specifically, the increase in $\mathcal{Q}_{\omega }$ can be attributed to three main factors:

1. (i) Boundary layers near tank walls: high velocity gradients close to the no-slip boundaries generate intense vorticity, contributing to enstrophy-related dissipation.

2. (ii) Fluid–solid impacts: as shown in previous sections, water impacts produce strong vorticity near the horizontal tank walls during the collapse of entrapped cavities, adding to $\mathcal{Q}_{\omega }$ .

3. (iii) Vortex structures from free-surface reconnection: the most significant contribution originates from vortices generated by the reconnection of fragmented free-surface elements, as discussed in § 5.4.3.

Figure 26.

Time history of the slosh dissipation energy with increasing spatial resolution for the $\delta$ -LES–SPH model. (a) without surface tension. (b) with surface tension.

Figure 26. Long description

Two line graphs depict the time history of slosh dissipation energy with increasing spatial resolution for the delta-LES-SPH model. Panel A shows the data without surface tension, while Panel B includes surface tension. Each graph features multiple lines representing different spatial resolutions, labeled as N equals 50, 100, 200, 400, and 800. The experimental data is also plotted for comparison. The x-axis in both graphs represents normalized time (t/T), and the y-axis represents normalized dissipation energy (Qdiss/Δε). In both panels, the dissipation energy decreases over time. Higher spatial resolutions generally show a closer fit to the experimental data. Without surface tension, the lines for higher resolutions (N equals 400 and 800) align more closely with the experimental data compared to lower resolutions. With surface tension, a similar trend is observed, but the lines for N equals 400 and 800 show a more pronounced decrease in dissipation energy over time, aligning well with the experimental data.

Figure 26 reports the non-dimensional total energy dissipation obtained from 2-D simulations at different spatial resolutions, both without and with surface tension (i.e. $\textit{We} = \infty$ and $\textit{We} = 2240$ ). For the case with surface tension, simulations were started at a minimum resolution of $N = 100$ , since lower resolutions were insufficient to capture the associated effects.

As emphasised in the previous section, the inherently chaotic nature of the flow required ten independent simulations for each spatial resolution to obtain statistically reliable results, both with and without surface tension. The only exception is the case with $N = 800$ , for which a single run was performed because of the prohibitive computational cost.

The results indicate that $N = 800$ practically represents the finest resolution attainable in the present study, with only limited changes expected under further refinement. For all resolutions up to $N \leq 400$ , the statistical dispersion associated with the independent realisations decreases as the number of particles increases. This indicates that the numerical solution is approaching a grid-converged 2-D limit.

The effect of surface tension ( $\textit{We} = 2240$ ) leads to a reduction in energy dissipation, associated with the formation of larger vortical structures and a smoother free surface. However, this stabilisation comes at the cost of a poorer agreement with the experimental data, suggesting that the influence of surface tension in a 2-D framework is significantly overestimated compared with the realistic 3-D dynamics.

Figure 27 presents the relative convergence rate of the energy-dissipation error for the 2-D simulations with and without surface tension. Since the experimental reference is inherently three-dimensional, these data cannot be used directly to compute the convergence rate. Nevertheless, the numerical results clearly show that, in the absence of surface tension, the convergence rate is approximately $1.2$ , while it increases to nearly $2$ when surface tension is included.

Figure 27.

Rate of convergence of the energy dissipation for the $\delta$ -LES–SPH model. (a) without surface tension. (b) with surface tension.

Figure 28.

Effect of spatial resolution on the energy dissipated by enstrophy, $\mathcal{Q}_\omega$ , for $\delta$ -LES–SPH simulations without (a) and with (b) surface tension.

Figure 28. Long description

Two line graphs compare energy dissipation in sloshing flows under different conditions. Panel A and Panel B each depict a line graph showing the normalized energy dissipation over time. The x-axis represents time normalized by period (t/T), and the y-axis represents the normalized energy dissipation (Q_omega/Delta_epsilon). Each graph includes multiple lines representing different values of N (50, 100, 200, 400, 800). Panel A shows data for We = infinity, while Panel B shows data for We = 2240. The lines indicate how energy dissipation changes over time for different spatial resolutions. In both panels, higher values of N generally show lower energy dissipation over time. The trends suggest that spatial resolution significantly affects energy dissipation in sloshing flows, with finer resolutions (higher N) leading to more accurate predictions of energy dissipation.

Figure 28 reports the cumulative dissipation of energy from enstrophy, $\mathcal{Q}_{\omega }$ , in 2-D simulations at different spatial resolutions. The left panel corresponds to the case without surface tension ( $\textit{We} = \infty$ ), while the right panel includes surface-tension effects ( $\textit{We} = 2240$ ). In both cases, enstrophy-related dissipation increases with spatial resolution, as finer particle sizes allow turbulent eddies to be captured more effectively.

However, the presence of surface tension reduces the overall dissipation of turbulent energy, since the free surface undergoes less fragmentation, leading to the generation of fewer vortical structures. Consequently, when surface tension is included, the enstrophy-related dissipation $\mathcal{Q}_{\omega }$ exhibits a clear convergence trend with increasing spatial resolution $N$ , with an estimated order of convergence close to unity. In contrast, in the absence of surface tension, $\mathcal{Q}_{\omega }$ does not show a convergent behaviour: each time $N$ is doubled, the dissipation increases significantly. This behaviour is associated with the fact that higher resolutions allow the formation of progressively smaller droplets during free-surface fragmentation, while simultaneously increasing their number. When these droplets impact the surrounding liquid – either underlying or overlying, depending on the stage of the tank motion – they generate additional small-scale vortices. To the authors’ knowledge, this behaviour has not been previously reported in the literature.

Figure 29 reports the time histories of the time-integrated enstrophy

(5.11) \begin{align} \mathcal{A}(t) := \int _{t_0}^{t} \int _{{D}} \| \boldsymbol{\nabla }\times \boldsymbol{u} \|^2 \, \text{d}V \, \text{d}t, \end{align}

normalised by the reference value $\varDelta \mathcal{E}/\mu$ , for the same 2-D simulations at five different spatial resolutions $N$ and for both cases, $\textit{We}=\infty$ and $\textit{We}=2240$ .

The first observation is that $\mathcal{A}(t)$ diverges almost linearly with the spatial resolution $N$ in both regimes, with the divergence being more pronounced in the absence of surface tension. This behaviour is physically expected and reflects the LES character of the present simulations: $\mathcal{A}(t)$ is dominated by the smallest resolved scales of $\| \boldsymbol{\nabla }\times \boldsymbol{u} \|^2$ , and at the highest resolution considered the in-plane discretisation $H/\Delta r = 800$ remains coarser than the Kolmogorov-scale ratio $H/\eta \simeq 2200$ (cf. § 5.4.1). Each further refinement therefore adds new contributions on previously unresolved scales, so that the integral grows without converging on the resolutions accessible in this work. As expected, the magnitude of $\mathcal{A}$ is reduced when surface tension is included, because of the less intense free-surface fragmentation and the resulting generation of fewer vortical structures.

Figure 29.

Effect of spatial resolution on the time-integrated enstrophy $\mathcal{A}$ for $\delta$ -LES–SPH simulations without (a) and with (b) surface tension.

It is essential to distinguish between $\mathcal{A}(t)$ , which has the dimensional structure of a small-scale-dominated quantity, and the physically meaningful dissipation $\mathcal{Q}_\omega$ . Considering the decomposition

(5.12) \begin{align} \mathcal{Q}_\omega = \mu \, \mathcal{A}(t) + \int _{t_0}^{t} \int _{{D}} \mu ^T \, \| \boldsymbol{\nabla }\times \boldsymbol{u} \|^2 \, \text{d}V \, \text{d}t, \end{align}

the eddy viscosity $\mu ^T \propto (C_S \Delta r)^2 \| \mathbb{D} \|$ provides the compensating factor that controls the divergence of $\mathcal{A}$ in the small-scale-dominated part of the integral: $\mu ^T$ decreases as $\Delta r$ decreases, in such a way that the SGS contribution to $\mathcal{Q}_\omega$ remains bounded. This is the expected behaviour of a well-posed LES closure and is precisely what produces the convergence of $\mathcal{Q}_\omega$ with surface tension reported above.

A quantitative measure of the role of the SGS closure is provided by the ratio of the resolved component $\mu \mathcal{A}$ to the total $\mathcal{Q}_\omega$ : at the highest resolution $N = 800$ , the resolved component is approximately $28\,\%$ of the total in both Weber-number regimes, so that the SGS model accounts for the remaining $\sim 72\,\%$ of the enstrophy-related dissipation. The same SGS-dominated regime applies to the total viscous dissipation: the resolved viscous dissipation $\mathcal{Q}_V$ accounts for approximately $22\,\%$ of the total $(\mathcal{Q}_V + \mathcal{Q}_V^{\textit{turb}})$ at $N = 800$ , in both Weber-number regimes. These ratios indicate that the simulations at $N = 800$ are well-resolved LES rather than DNS, and that the SGS model plays a quantitatively important role in the dissipation budget. The distinction between $\mathcal{A}(t)$ , which grows with $N$ as expected in LES, and $\mathcal{Q}_\omega$ , which converges with surface tension once the SGS contribution is included, is therefore not a contradiction but a direct manifestation of the LES closure at work.

Figure 30 shows the cumulative energy dissipation with respect to different components at the highest resolution $N=800$ , highlighting the effect of surface tension ( $\textit{We}=2240$ , right panel) in limiting free-surface fragmentation and modifying the distribution of dissipative events. The left panel corresponds to the case without surface tension ( $\textit{We}=\infty$ ), where droplets progressively decrease in size because of the absence of a physical constraint, resulting in higher total energy dissipation and increased generation of small-scale vortices. The right panel illustrates the case with surface tension, where larger droplets and more coherent jets reduce overall energy dissipation, limit small-scale vortices and decrease deformation of the fluid domain.

Comparing the two cases, $\mathcal{Q}_\omega$ , expressed as a percentage, remains approximately $30\,\%$ despite the lower total energy dissipation with surface tension, while $\mathcal{Q}_N$ , associated with droplet impacts, is higher with surface tension because of stronger impacts from larger droplets and more coherent jets. Moreover, $\mathcal{E}_{\textit{Def}}$ is lower with surface tension, reflecting reduced deformation and smoother evolution of the free surface. These results demonstrate that surface tension significantly modulates the distribution of dissipative mechanisms, reducing fragmentation while altering impact dynamics and enstrophy-related dissipation.

Figure 30.

Decomposition of the slosh dissipation for the highest resolution $N=800$ (see (3.6)) without (a) and with (b) surface tension.

Figure 30. Long description

Two line graphs depict the decomposition of slosh dissipation for the highest resolution. Panel A on the left shows the data without surface tension, while Panel B on the right includes surface tension. Each graph features a vertical axis labeled Q/Δε and a horizontal axis labeled t/T. The graphs include multiple lines representing different components: Qω, (QV + Qturb), Qdiss, and QN. In Panel A, Qω is represented by an orange line, (QV + Qturb) by a blue line, and Qdiss by a dashed gray line. The graph shows that Qω remains relatively constant, while (QV + Qturb) and Qdiss decrease over time. Key percentages are marked, such as 45 percent for εDef and 25 percent for QN. Panel B follows a similar structure with Qω, (QV + Qturb), and Qdiss lines. Here, Qω and (QV + Qturb) also remain relatively constant, while Qdiss decreases. Notable percentages include 29 percent for both Qω and εDef, and 42 percent for QN. The graphs illustrate the impact of surface tension on the dissipation dynamics over time.

These observations confirm that, in the present violent-sloshing regime, the dissipation budget is shaped primarily by fragmentation and reconnection of the free surface, with a secondary modulation introduced by surface tension. The persistence of a nearly constant enstrophy share ( $\mathcal{Q}_\omega \approx 30\,\%$ of the total) across the two fragmentation regimes considered here is an empirical statistical signature of this distribution, rather than a claim about a causal mechanism: in the present class of flows, the generation and destruction of coherent vortical structures is dominated by interface-driven events, not by vortex stretching as in canonical 3-D turbulence. The redistribution between droplet-impact dissipation ( $\mathcal{Q}_N$ ) and deformation-related dissipation ( $\mathcal{E}_{\textit{Def}}$ ) reflects the structural changes induced by surface tension.

From a numerical perspective, the 2-D approach captures not only the qualitative trends but also the global integral quantities, such as the vertical force on the tank and the total dissipated energy, within the dispersion of the chaotic flow. This is consistent with the close agreement between 2-D and 3-D predictions documented in the literature for this class of bounded violent free-surface flows (Lubin et al. Reference Lubin, Vincent, Abadie and Caltagirone2006; Marrone et al. Reference Marrone, Colagrossi, Gambioli and González-Gutiérrez2021b ) and directly verified through the reduced-width 3-D simulations of § 5.3. The local statistics associated with the lifespan of coherent vortical structures are likewise consistent between the 2-D and 3-D descriptions. The 2-D analysis is therefore not an approximation of the 3-D case in this regime, but a faithful description of the dynamics that controls the engineering behaviour of the system.

In summary, figure 31 reports the total energy dissipation for both the 2-D ( $N=800$ ) and 3-D ( $N=75$ ) cases, compared against the experimental data. In this plot, surface tension is neglected (i.e. $\textit{We}=\infty$ ): as already emphasised, in the 3-D framework the spatial resolution is insufficient to capture surface-tension effects, whereas in the 2-D setting such effects tend to be overestimated. Notably, the comparison shows that both 2-D and 3-D simulations, despite similar total particle counts (approximately $2.2\times 10^6$ in two dimensions and $2.8\times 10^6$ in three dimensions), yield numerical results close to the experimental measurements. A plausible explanation is linked to the mechanisms discussed above: in violent sloshing, vortex generation, merging and break-up are strongly constrained by free-surface fragmentation and the reconnection dynamics.

As a result, the effective energy-transfer pathways in two and three dimensions become less dissimilar than in canonical turbulence, which helps rationalise the close agreement in the total dissipated energy obtained with comparable particle counts. Overall, the present results indicate that, in violent sloshing, the free surface is an active ingredient of the turbulence dynamics, as it both generates vorticity through impacts and reconnection and limits coherent-eddy lifetimes through interface-induced break-up.

Figure 31.

Comparison against the experimental data of the slosh dissipation, $\mathcal{Q}_{\textit{diss}}$ , evaluated through $\delta$ -LES–SPH simulations in the 2-D and 3-D frameworks. Surface tension is not considered, i.e. $\textit{We}=\infty$ .