Numerical Calculation of Slosh Dissipation

Abstract

As part of the Sloshing Wing Dynamics H2020 EU project, an experimental campaign was conducted to study slosh-induced damping in a vertically excited tank filled with liquid water or oil and air. In this work, we simulate these experiments using two numerical approaches. First, a single-phase, weakly compressible liquid model is used, and the gas flow (air) is not modeled. For this approach, a proven Smoothed Particle Hydrodynamics (SPH) model is used. In the second approach, both phases are simulated with an incompressible liquid and weakly compressible gas model via a Finite Volume Method (FVM) using Volume-of-Fluid (VOF) to track the liquid phase. In both approaches, the energy distribution of the flow is calculated over time in two- and three-dimensional simulations. It is found that there is reasonable agreement on the energy dissipation evolution between the methods. Both approaches show converging results in 2D simulations, although the SPH simulations seem to have a faster convergence rate. In general, the SPH results tend to overpredict the total dissipation compared to the experiment, while the finite volume 2D results underpredict it. Time histories of the center of mass positions are also compared. The SPH results show a much larger vertical center of mass motion compared to the FVM results, which is more pronounced for the high Reynolds number (water) case, probably linked to the absence of the air phase. On the other hand, the limited center of mass motion of the FVM could be linked to the need for higher spatial resolutions in order to resolve the complex gas–liquid interactions, particularly in 3D.

1. Introduction

The Sloshing Wing Dynamics Project (https://slowd-project.eu/, accessed on 16 September 2019) aims to study fluid sloshing in commercial aircraft wing tanks and how wing loads can be reduced as a result [1].

The project includes a future experimental campaign of a full scale aircraft wing structure and corresponding internal fuel tanks. The intention is to excite the structure with a variety of loading conditions in which slosh loads are the largest.

Various experimental studies were performed during the project [2,3,4] on a smaller scale. Some of the findings from these studies included that the two-phase flow in the tank contributes significantly to damping wing vibrations. The flow in wing tanks is particularly violent in many load cases of interest. The liquid and gas phases are mixed well and contain small structures (films, ligaments, drops and bubbles) with the flow often undergoing a significant amount of liquid–wall and liquid–liquid impact events during the tank excitation.

In spite of its industrial impact, there are relatively few studies dealing with vertically excited sloshing flows in the literature. Experimentally, one of the first studies dealing with this problem is [5] and, more recently [6]. A model scale experimental campaign for aircraft wing excitation is carried out in [7]. In these experiments, a partially filled tank is fixed rigidly to an elastic beam, which is deflected to induce a heaving oscillatory motion with accelerations up to 10 g. In [8,9], experiments of slosh-induced energy dissipation in partially filled tanks subject to vertical, harmonic excitation is investigated. In [10] the decaying harmonic excitation is investigated experimentally, along with an equivalent mechanical model and a weakly compressible single-phase SPH simulations.

In the present work, we perform numerical simulations of the single degree of freedom vertical slosh experiments carried out in [4]. This experimental work has recently been extended [11] with comments on repeatability as well as a study on scaling effects. Previous results on the experiments in [4] are published [12,13], but new results are added in this work. Similar to [12,13], we employ an energetic approach to study the dissipation mechanisms in the tank flow. This is completed for two different numerical approaches. In the first approach, only the liquid is solved, using a weakly compressible formulation with a Smoothed Particle Hydrodynamics (SPH) method. The second approach uses the multi-physics code Elemental^® [14,15,16,17,18,19], which employs a vertex-centered, Finite Volume Method (FVM) with Volume-of-Fluid (VOF) tracking of the liquid phase. Both the liquid and gas (air) are therefore solved. The liquid is assumed incompressible and the gas weakly compressible [20]. Surface tension effects are accounted for using a Continuum Surface Force (CSF) method, where the SPH approach does not include capillary effects. As explained in [21,22], the latter does not have a significant role in the fluid energy dissipation mechanisms (apart from the initial liquid meniscus). This is due to the high Weber numbers of the simulated flows. As noted, surface tension is indeed relevant during the initial upward acceleration of the tank when a Rayleigh–Taylor instability develops (see, e.g., [23]). The instability is influenced by the fluid meniscus formed at the intersection between the lateral walls and the liquid surface as the accelerating meniscus perturbs the interface, leading to the growth of a Rayleigh–Taylor instability during the upward acceleration of the tank from rest. For this reason, the initial meniscus free-surface deformation is taken into account in the simulations (see also, e.g., [12,22]).

Test cases are set up and solved using similar initial conditions and tank acceleration profiles for both approaches. The energy distribution in both cases is studied, as well as other metrics, such as the center of mass position. In this work, we prescribe the tank acceleration to match experimental measurements. This allows us a more direct comparison to the experiment and ensures that both methods are subject to the same acceleration. There has been research with coupled simulation, where the fluid–structure interaction is resolved [13], but that is outside the scope of this work. First, two-dimensional simulations are performed using random interface perturbations to study the variation in energy metrics with the slight difference in interface shape. Three-dimensional experiments are also performed to compare the 2D results of each approach as well as an overall comparison between approaches. The results of these cases provide insights into several aspects. These include the suitability of the modeling approach (single or two-phase), the chaotic nature of the flow (in 2D) with random initial perturbations, mesh resolution effects, and the difference in 3D results compared to 2D.

The experiment will be described next, followed by the mathematical formulation and numerical method details of each approach.

Problem Description

The experimental campaign on which simulations are based is detailed in [4]. A tank is connected to a set of six springs: three on the upper and three on the lower side. The rectangular tank has the following size: $L=10$ cm (length), $D=6$ cm (height), $W=6$ cm (depth), while the filling height is denoted H. The tank is filled up to $50\%$ of its volume (i.e., $H=3$ cm), and water and oil are considered as sloshing liquids (water mass of $m_l$ = $0.18$ kg and oil $m_l=0.162$ kg). This is to study both low and high Reynolds number cases.

When the springs are released, the tank oscillates at a characteristic frequency of about $f_0=1/T$ = 6.51 Hz for both the liquids. Figure 1 sketches three pictures taken at different time instants during the experimental test (more details on the experiment can be found in the article by [4]). At time $t=T$, the tank reaches its maximum vertical displacement $2A=1.14L$ and a maximum acceleration $a_{max}\approx$ 8.8 g. The latter value is close to 10 g, which is of the same order of magnitude as limit design cases of an airliner wing during a wind gust (see [2]). The maximum velocity $U_{max}$ of about 2.33 m/s is reached by the tank at $t=0.747T$ during the first rising stage. The reference energy adopted is the potential energy (of all fluid in the tank in its initial configuration) when the tank reaches its maximum vertical displacement, i.e., $\Delta E=m_{l}g2A$, which is equal to $\Delta E_w=0.2013$ J for the water and $\Delta E_o=0.1812$ J for the oil. The reference power is $\Delta P=\Delta E/T$.

The Reynolds number of the problem is calculated using the maximum tank velocity and length dimension (Re = $U_{\max }L/\nu$), where $\nu$ is the kinematic viscosity of the liquid. The resulting Reynolds numbers are 233,000 for water and 4660 for oil.

Figure 2 shows the recorded vertical motion of the tank considering both water and oil. In the oil test case, the total amount of dissipated energy is lower compared to water. This counter-intuitive result is, indeed, in agreement with the experimental findings in [24] for a tuned liquid damping system. It is explained that water sloshing flows dissipate energy mostly through energetic breaking wave phenomena. In liquids characterized by higher viscosity, the breaking events are less intense and, therefore, less energy can be dissipated. An analogous result was found in [24,25], where the authors underline that the higher fragmentation phenomena occurring at higher Reynolds number induces a larger fluid energy dissipation. This was similarly found in [22], where greater fragmentation resulted in greater liquid impact.

An overview of the respective methods, their mathematical formulations and numerical methods will be provided next, starting with the SPH method.

2. SPH Method: The SPH-Flow Solver

2.1. Governing Equations Adopted for the SPH Model

The main simplifications adopted for the SPH model are:

• Modeling of only the liquid phase;
• Neglect of thermal conductivity and surface tension;
• Liquid is assumed to be a weakly compressible media through an artificial speed of sound.

The use of a single-phase model obliges the neglect of all physics related to the air contained in the fluid domain. This assumption may appear inappropriate for violent sloshing simulations where air entrapment is unavoidable. However, in [26], it was shown that although the air phase plays a relevant role in the flow evolution, the evaluation of the energy dissipation in violent flows, even under the single-phase hypothesis, is still accurate enough. A further validation can also be found in [24,25], where the study of the mechanical energy dissipation induced by sloshing and wave breaking in a fully coupled angular motion system was investigated. In this work, a single-phase SPH model is demonstrated to be able to predict the experimental results of the fluid dissipated energy with reasonable accuracy.

The role of surface tension is negligible according to the velocity and length scales considered, as also explained in [21], since the Weber number, $\mathrm{We}=\rho U_{\max }^{2}D/\sigma$ (where $\sigma$ is the surface tension coefficient of the liquid considered) is larger than 4000 for both the liquids considered in this work. However, the liquid menisci due to the contact angle between the lateral walls and the liquid surface need to be considered in the initial conditions of the simulations presented in this work, since they form the main perturbation mechanism of the interface leading to the Rayleigh–Taylor instability during the initial flow stage (see also [23]).

In the SPH model, the flow evolution is considered governed by the weakly compressible Navier–Stokes equations:

$$\left\{\begin{array}{cc}{\displaystyle \frac{D\rho }{Dt}=-\rho \mathrm{d}\mathrm{i}\mathrm{v}\left(\mathit{u}\right),}\hfill & {\displaystyle \frac{D\mathit{u}}{Dt}=g+f_{NI}+\frac{\mathrm{d}\mathrm{i}\mathrm{v}\left(\tau \right)}{\rho }}\hfill \\ {\displaystyle \frac{De}{Dt}=\frac{:D}{\rho },}\hfill & {\displaystyle \frac{D\mathit{r}}{Dt}=\mathit{u},p=c_0^2\left(\rho -{\rho }_0\right),}\hfill \end{array}\right.$$

(1)

where $D/Dt$ represents the Lagrangian derivative, u the flow velocity, r the position of the material points, $\rho$ the fluid density, $\tau$ the stress tensor, D the rate of strain tensor and $\mathit{g}$ the gravitational acceleration. Regarding the equation of state, which links pressure to the density field, $c_0$ is the speed of sound of the liquid, while the reference density value ${\rho }_0$ refers to the density along the free surface (where p is assumed to be equal to zero). Under the weakly-compressible hypothesis, the effects of entropy/temperature on the pressure are assumed to be negligible. Therefore, the pressure p is assumed to be linearly dependent on the density and proportional to the speed of sound of the liquid $c_0$. The weakly-compressible condition requires that:

$$c_0\gg max\left(U_{\max },\sqrt{\frac{{(\Delta p)}_{\max }}{\rho }}\right),$$

(2)

where $U_{\max }$ and ${(\Delta p)}_{\max }$, respectively, are the maximum fluid speed and the maximum pressure variation expected (with respect to the zero pressure free-surface level) within the fluid domain. To avoid too small time steps, $c_0$ is always set lower than its physical counterpart (in the present work, about two orders of magnitude lower). The constraint (2), however, must be respected in the interest of the validity of the weakly compressible assumption. As discussed in [27,28], if this weakly compressible constraint is indeed satisfied, the energy dissipation linked to water impacts using Equation (1) is consistent with that of incompressible flow models.

The thermal conductivity and the surface tension effects are neglected, and the liquid is assumed to be Newtonian, i.e.,: $\tau =[-p+\lambda \mathrm{d}\mathrm{i}\mathrm{v}(u\left)\right]I+2\mu D$, where $\mu$ and $\lambda$ are the primary and secondary dynamic viscosities of the liquid, and I the identity tensor.

Equation (1) can be either solved in an Inertial Frame of Reference (I-FoR) where the tank is moving (relative to the reference frame) or in a Non-inertial Frame of Reference (Ni-FoR), in which the reference frame is fixed to the tank. The Ni-FoR, adopted in this work, is generally preferred due to its easier implementation. It uses a static mesh with tank accelerations applied as an acceleration source term $f_{NI}$ to the fluid in the governing Equation (1).

2.2. Brief Recall of the $\delta$-LES-SPH Scheme

In the present section, we briefly introduce the SPH scheme defined in [29,30]. The governing Equation (1) is discretized as a set of fluid particles whose masses $m_i$ remain constant during the motion. The particles are set initially on a lattice with homogeneous spacing $\Delta r$, and hence, the particles’ volumes $V_i$ are evaluated initially as $\Delta r^n$ (where n is the number of spatial dimensions). The particle masses $m_i$ are calculated through the initial density field using the equation of state and the initial pressure field. During the time evolution, volumes $V_i$ change in time according to the particle density ${\rho }_i$.

For the sake of brevity, in the following text the notation $r_{ji}$ indicates the differences in the particles positions $(r_j-r_i)$, and the same holds for the velocity fields $u_{ji}$ and $\delta u_{ji}$, while, for the generic scalar field the notation $f_{ij}$ just indicates the dependency of the field f on the indices i and j. The spatial gradients are approximated through convolution summations with a kernel function $W_{ij}$. This function has a compact support whose reference length is denoted by h and referred to as the smoothing length. As in [29], a C2-Wendland kernel is adopted in the present work. For this kernel, the radius of the support is $2h$. In the numerical simulations that follow, the number of neighboring particles is chosen as $2h=4\Delta r$ in two dimensions and $2h=3\Delta r$ in three dimensions.

The governing Equation (1) is discretized within the $\delta$-LES-SPH scheme as:

$$\left\{{\displaystyle \begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}}& {\displaystyle =-{\rho }_i\sum _j(u_{ji}+\delta u_{ji})·{\nabla }_i{W}_{ij}{V}_j+\sum _j({\rho }_j\delta u_j+{\rho }_i\delta u_i)·{\nabla }_i{W}_{ij}{V}_j+{\mathcal{D}}_i^{\rho }}\hfill \\ \hfill {\rho }_i\frac{\mathrm{d}{u}_i}{\mathrm{d}t}& =F_i^p+F_i^v+{\rho }_i(g+f_{NI})+{\displaystyle \sum _j}({\rho }_j{u}_j\otimes \delta u_j+{\rho }_i{u}_i\otimes \delta u_i)·{\nabla }_i{W}_{ij}{V}_j\hfill \\ \frac{\mathrm{d}{r}_i}{\mathrm{d}t}& =u_i+\delta u_i,V_i=m_i/{\rho }_i,p=c_0^2(\rho -{\rho }_0),\hfill \end{array}}\right.$$

(3)

where the indexes i and j refer to generic i-th and j-th particles, $F_i^p$ and $F_i^v$ are the pressure and viscous forces acting on the particle i. The vector $\delta u$ is the particle shifting velocity adopted to regularize the particles’ spatial distribution during their motion. The time derivative $d/dt$ used in (3) indicates a quasi-Lagrangian derivative since the particles are moving with the modified velocity $(u+\delta u)$, and the above equations are written in an Arbitrary Lagrangian–Eulerian framework. For this reason, the continuity and the momentum equations contain terms with spatial derivatives of $\delta u$ (for details, see [31]).

The term ${\mathcal{D}}_i^{\rho }$ is the numerical diffusive term introduced by [32] to filter out the spurious high-frequency noise in the pressure field:

$${\mathcal{D}}_i^{\rho }:=\sum _j{\delta }_{ij}{\psi }_{ji}·{\nabla }_i{W}_{ij}{V}_j$$

(4)

The law adopted in the present work for the particle shifting velocity $\delta u$ term, as well as the functions ${\delta }_{ij}$ and ${\psi }_{ij}$, are given in [12] and not reported here for the sake of brevity.

Regarding the pressure and the viscous forces, these are expressed as:

$$\left\{\begin{array}{ccc}{F}_i^p\hfill & :=\hfill & {\displaystyle -\sum _j\left(p_j+p_i\right){\nabla }_i{W}_{ij}{V}_j}\hfill \\ {\displaystyle F_i^v}\hfill & :=\hfill & {\displaystyle K\sum _j(\mu +{\mu }_{ij}^T){\pi }_{ij}{\nabla }_i{W}_{ij}{V}_j}\hfill & \hfill & {\displaystyle K:=2(n+2)}\hfill \\ {\displaystyle {\pi }_{ij}}\hfill & :=\hfill & {\displaystyle \frac{u_{ij}·r_{ij}}{\left|\right|r_{ji}{\left|\right|}^2}{\mu }_{ij}^T:=2\frac{{\mu }_i^T{\mu }_j^T}{{\mu }_i^T+{\mu }_j^T}}\hfill & \hfill & {\displaystyle {\mu }_i^T:={\rho }_0{\left(C_{S}l\right)}^2\left|\right|D_i\left|\right|}\hfill \end{array}\right.$$

(5)

where $C_S$ is the so-called Smagorinsky constant set equal to 0.18 (see [33]). The viscous term (5) contains both the effect of the real viscosity $\mu$, as well as the one related to the local turbulent viscosity ${\mu }_i^T$.

$\delta$-LES-SPH Scheme: Enforcement of the Boundary Conditions

The kinematic and dynamic free-surface and no-slip conditions are applied, respectively, on the free surface and on the solid surfaces and thus coupled with the governing Equation (1). More specifically, while the free surface boundary conditions are intrinsically satisfied in the SPH methods (see [34]), it is difficult to resolve thin wall boundary layers (WBL) as easily performed in mesh-based methods. Concerning the latter, if a solution is represented by the use of smaller particles close to the walls, it implies large CPU costs linked to the explicit time integration of the scheme. Indeed, in high Reynolds number conditions, thin WBLs are developed and a simple no-penetration boundary condition (free-slip) is preferred to avoid a too-demanding cost in terms of computational resources. The Reynolds number of the simulations with water, presented in this work, is about 233,000 and an initial estimation of the WBL thickness returns in that the near-wall regions, even with the maximum spatial resolution, are still under-resolved so the free-slip conditions are preferred. On the other hand, the Reynolds number in the simulations with oil is about 4660, and in this condition, the particle size used for water is sufficient to resolve the boundary layer developed by the oil, thus allowing the application of the no-slip conditions. The limit of the free-slip assumption for the water simulations is evaluated in [12], and it returns that the energy dissipation is not substantially influenced by choice of the wall boundary conditions in the violent free-surface flows presented here. This result is also connected to two opposing phenomena: the wall friction slows down the run-up of the liquid jet, decreasing the intensity of the liquid impact against the tank walls, but it represents a further energy dissipation mechanism.

2.3. $\delta$-LES-SPH Scheme: Evaluation of the Energy Dissipation

Following the analysis performed in [30,35], the energy balance for the particle system is extended to the $\delta$-LES-SPH equations. For the sake of brevity, only the main terms are briefly reported in this section. The $\delta$-LES-SPH energy balance can be written as:

$$\left\{\begin{array}{c}{\displaystyle {\dot{E}}_K+{\dot{E}}_P-{\mathcal{P}}_{NF}={\mathcal{P}}_V+{\mathcal{P}}_V^{turb}+{\mathcal{P}}_{\mathcal{N}}}\hfill \\ {\displaystyle E_K\left(t\right)=\frac{1}{2}\sum _i{m}_i{u}_i^2,E_P\left(t\right)=\sum _i{m}_{i}g{y}_i,{\mathcal{P}}_{NF}=\sum _i{m}_i(f_{NI}·u_i)}\hfill \end{array}\right.$$

(6)

where on the left-hand side, $E_K$ and $E_P$ are the kinetic and potential energy of the particle system, and ${\mathcal{P}}_{NF}$ is the power linked to the non-inertial forces. In Equation (6), $y_i$ is the vertical position of the generic i-th particle.

Conversely, the right-hand side of the energy balance (6) contains the dissipation terms due to the real viscosity ${\mathcal{P}}_V$, to the turbulent viscosity ${\mathcal{P}}_V^{turb}$, while ${\mathcal{P}}_N$ takes into account the effect of the density diffusion and the particle shifting $\delta u$.

The elastic potential energy linked to the compressibility of the liquid is negligible within the weakly-compressible assumption; hence, it is not considered in the energy balance.

The power related to the viscous forces is directly evaluated through the expressions (5) as:

$${\mathcal{P}}_V+{\mathcal{P}}_V^{turb}=\frac{K}{2}\sum _i\sum _j(\mu +{\mu }_{ij}^T){\pi }_{ij}{u}_{ij}·{\nabla }_i{W}_{ij}{V}_i{V}_j$$

(7)

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

Finally, the energy dissipated by the fluid can be evaluated by integrating time into Equation (6)

$$\left\{\begin{array}{c}{\displaystyle E_K\left(t\right)+E_P\left(t\right)-W_{NF}\left(t\right)=E_{diss}\left(t\right),}\hfill \\ {\displaystyle W_{NF}\left(t\right)={\int }_{t_0}^t\left(\sum _i{m}_i{f}_{NI}·u_i\right)dt,E_{diss}\left(t\right)={\int }_{t_0}^t\left({\mathcal{P}}_V+{\mathcal{P}}_V^{turb}+{\mathcal{P}}_{\mathcal{N}}\right)dt}\hfill \end{array}\right.$$

(8)

where $W_{NF}$ is the work performed by the non-inertial forces on the fluid. The first equation of (8) gives two ways for the evaluation of the energy dissipated by the fluid: (i) the first is using the left-hand side, while (ii) a second way is to use the second equation in which the three dissipation terms are integrated in time. Both these approaches were adopted in the SPH simulations in order to verify that the present model is able to close the energy balance accurately.

3. Relation between Energy Dissipation and Fluid Center of Mass Motion

Considering the motion of the tank as a pure vertical translation, following the analysis in [12], the energy balance can be reshaped as:

$$E_K\left(t\right)+m_{l}g{y}_G\left(t\right)+m_l{\int }_{t_0}^t{\dot{y}}_G\left(t\right)a_{tank}\left(t\right)dt=E_{diss}\left(t\right),y_G\left(t\right):=\frac{{\sum }_i{\rho }_i{V}_i{y}_i}{m_l},$$

(9)

where $a_{tank}$ is the vertical acceleration of the tank and $y_G$ the position of the fluid center of mass in the Ni-FoR. The latter is evaluated in a discrete form on the finite volume cells or on the SPH particles.

The first equation of (9) highlights the role of the motion of the fluid center of mass in the Ni-FoR concerning the energy dissipated by the liquid sloshing. In particular, the phase lag between the velocity ${\dot{y}}_G$ and $a_{tank}\left(t\right)$ is a crucial aspect of the phenomenon. The third term of the left-hand side of the first Equation of (9) is the work completed by the non-inertial force acting on the liquid:

$$W_{NF}\left(t\right)={\int }_{t_0}^t{\mathcal{P}}_{NF}dt=-m_l{\int }_{t_0}^t{\dot{y}}_G\left(t\right)a_{tank}\left(t\right)dt$$

(10)

In the problem studied in this article, $W_{NF}$ is much larger than the kinetic and potential energy in the Ni-FoR, and therefore, as the time increases, $W_{NF}$ becomes increasingly closer to the dissipated energy $E_{diss}$. Considering the energy balance (9) at the time $t_f=30$ T at the end of the experiment, when the fluid almost comes back to a rest condition, they are going to be very close to each other:

$$W_{NF}\left(t_f\right)\approx E_{diss}\left(t_f\right)$$

Following again [12], one can find that the acceleration of the center of mass ${\ddot{y}}_G$ is linked to the vertical forces $F_y$ that the liquid exerts on the tank:

$$m_l{\ddot{y}}_G\left(t\right)=-F_y\left(t\right)-m_l[g+a_{tank}\left(t\right)]$$

(11)

and the work completed by those forces of the above expression is:

$$W_{ext}^{dyn}\left(t\right)=m_l{\int }_{t_0}^t{\ddot{y}}_G\left(t\right)v_{tank}\left(t\right)dt$$

(12)

If the fluid is “frozen” (i.e., no liquid sloshing), the above external work performed by the tank walls on the fluid is zero since the center of mass $y_G$ is a constant in the Ni-FoR. Considering that at $t=t_0$, the tank is in a rest condition, using integration by parts, it follows that:

$$W_{ext}^{dyn}\left(t\right)=m_l{\dot{y}}_G\left(t\right)v_{tank}\left(t\right)+W_{NF}\left(t\right)$$

(13)

where, as explained in [12], the first term on the right-hand side is linked to differences in the mechanical energy (i.e., $E_K+E_P$) evaluated in the I-FoR and in the Ni-FoR. Again, similarly to what was stated above, at the end of the experiment (i.e., $t=t_f$), the tank does not move anymore; hence, the two work terms $W_{ext}^{dyn}$ and $W_{NF}$, coincide and are practically both equal to the dissipated energy $E_{diss}$ (see also Section 5.1).

It is worth noting that the center of mass motion $y_G\left(t\right)$ can be experimentally evaluated measuring $F_y$ with a load cell, while the recording of $y_G\left(t\right)$ from video-camera frames would require a more complex indirect measure, especially in a three-dimensional context [private communication, Gambioli, 2021] [36]. Finally, all the relations written in this section are valid both for no-slip and free-slip conditions on the tank walls. In fact, those boundary conditions influence the fluid motion, which reflects directly on the center of mass motion $y_G\left(t\right)$.

4. Volume-of-Fluid Finite Volume Method: The Elemental Solver

4.1. Background

Consider the generic domain V, which is undergoing a vertical, rigid body acceleration ${\mathit{a}}_d$ and contains a liquid–gas mixture. A schematic illustration of the inertial reference frame $X,Y,Z$ and domain volume V inside non-inertial reference frame $x,y,z$ is shown in Figure 3.

The absolute velocity ${\mathit{u}}_a$ of the fluid in the global, inertial reference frame is defined as

$${\mathit{u}}_a(X,Y,Z,t)=\mathit{u}(x,y,z,t)+{\mathit{u}}_d(X,Y,Z,t),$$

(14)

where ${\mathit{u}}_d$ is the velocity of the fluid domain relative to the inertial reference frame ($X,Y,Z$) and $\mathit{u}$ is the fluid velocity relative to the domain (non-inertial reference frame and simulation domain $x,y,z$).

4.2. Governing Equations for the Elemental Vertex-Centered FVM

In this section, we present the governing equations used for the Elemental^® vertex-centered, finite volume VOF solver [19]. They consist of equations for the conservation of mass, momentum and an advection equation for the VOF, $\alpha$. The equations are written in integral form for a volume $V_i$, which represents a finite volume cell in tank volume V

$$\begin{array}{cc}\hfill \frac{1}{V_i}{\int }_{V_i}\left[\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·\rho \mathit{u}\otimes \mathbf{u}\right]\mathrm{d}V& =\frac{1}{V_i}{\int }_{V_i}\left[\nabla ·\mathit{\tau }+\rho \left(\mathit{g}+\mathit{a}\right)+\sigma \kappa {\mathit{n}}_{\alpha }{\delta }_{\alpha }\right]\mathrm{d}V\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{1}{V_i}{\int }_{V_i}\nabla ·\mathit{u}\mathrm{d}V& =\frac{1}{V_i}{\int }_{V_i}\frac{(\alpha -1)}{{\rho }_g}\frac{1}{c_g^2}\frac{\partial p}{\partial t}\mathrm{d}V,\hfill \end{array}$$

(16)

with $\rho$ the density, $\mathit{u}$ the fluid velocity relative to the tank, $\mathit{g}$ gravitational acceleration and $\mathit{a}$ the acceleration of the tank. The last term in (15) is the estimation for the surface force due to surface tension, with ${\delta }_{\alpha }$ a Dirac delta function defined on the VOF interface, ${\mathit{n}}_{\alpha }$ the interface normal, $\sigma$ the surface tension coefficient and $\kappa$ the interface curvature. In (16), a weakly compressible gas formulation [20] is used, with $\alpha$ the VOF, which is defined in volume $V_i$ as a volume weighted tracer of the liquid phase $H_{\alpha }$, so that $H_{\alpha }=1$ is in the liquid and $H_{\alpha }=0$ is in the gas phase

$${\alpha }_i=\frac{1}{V_i}{\int }_{V_i}{H}_{\alpha }\mathrm{d}V.$$

(17)

Further, $c_g$ is the acoustic velocity in the gas phase. The fluid density and viscosity are represented as $\rho$ and $\mu$. Throughout this work, a harmonic mean is used for the fluid viscosity, assuming constant values in each phase. The density within a finite volume $V_i$ is calculated as an average using the VOF

$$\rho =\alpha {\rho }_l+(1-\alpha ){\rho }_g,$$

(18)

where the liquid density ${\rho }_l$ is a constant in this work, and the gas density ${\rho }_g$ is calculated assuming an ideal gas law with adiabatic compression such that

$${\rho }_{g_2}={\rho }_{g_1}{\left(\frac{p_2}{p_1}\right)}^{(1/\gamma )}.$$

(19)

Here, $\gamma =\frac{c_p}{c_v}=1.4$ for air, with $c_p$ and $c_v$, respectively, the specific heat capacity at constant pressure and volume. The gas acoustic velocity is

$$c^2=\frac{\gamma p}{{\rho }_g},$$

(20)

where for dry air at Normal Temperature and Pressure (NTP) $\rho$ = $1.205$ kg·m${}^{-3}$ and $p=101.325$ kPa, hence $c_g=343$ m·s${}^{-1}$.

4.3. Numerical Method

Equations (15) and (16) are solved using a fractional step projection method [19]. First, a predicted velocity ${\mathit{u}}^{*}$ is calculated from (15), without the pressure term from the values at time step n:

$$\begin{array}{cc}{\displaystyle {\mathit{u}}^{*}={\displaystyle \frac{1}{{\rho }^{n+1}}}[}\hfill & {\displaystyle {\rho }^n{\mathit{u}}^n+\Delta t(-\nabla ·\rho \mathit{u}\otimes \mathit{u}{{|}^n+\nabla ·(\mu +{\mu }_t)\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)|}^n+}\hfill \\ \hfill & {\displaystyle {+\rho \mathit{g}|}^n{+\sigma \kappa \nabla \alpha |}^n\left)\right].}\hfill \end{array}$$

(21)

where ${\mu }_t$ is the turbulent viscosity which is computed using the LES Smagorinsky–Lilly model (as per the SPH method).

The pressure term contribution is added to obtain the velocity at the following time step $n+1$:

$${\mathit{u}}^{n+1}={\mathit{u}}^{*}-{\displaystyle \frac{\Delta t}{{\rho }^{n+1}}}\nabla p^{n+1}$$

(22)

Mass conservation (16) is then applied to this equation so that

$$\nabla ·{\mathit{u}}^{n+1}=\nabla ·\left({\mathit{u}}^{*}-{\displaystyle \frac{\Delta t}{{\rho }^{n+1}}}\nabla p^{n+1}\right)=-\frac{\left(1-{\alpha }^{n+1}\right)}{{\rho }_g{c}_g^2}{\displaystyle \frac{p^{n+1}-p^n}{\Delta t}}$$

$$\therefore 0=\nabla ·\left({\mathit{u}}^{*}-{\displaystyle \frac{\Delta t}{{\rho }^{n+1}}}\nabla p^{n+1}\right)+\frac{\left(1-{\alpha }^{n+1}\right)}{\gamma p^{n+1}}{\displaystyle \frac{p^{n+1}-p^n}{\Delta t}}.$$

(23)

The VOF is updated in each cell using a compressive, algebraic advection scheme [37], which ensures mass conservation [38]:

$${\alpha }^{n+1}={\alpha }^n-\Delta t\left(\nabla ·{\mathit{u}}^n{\alpha }^n\right)$$

(24)

A notionally second-order predictor-corrector time integration scheme was used in this work, so that Equations (21)–(24) are iterated twice per time step.

4.4. Mechanical Energy Calculation: Elemental

The momentum conservation of the fluid system can be written in the inertial reference frame as:

$$\frac{D}{Dt}\left(\rho {\mathit{u}}_a\right)=\nabla ·\mathbf{\tau }+\rho \mathit{g}+\sigma \kappa \nabla \alpha$$

(25)

where $D/Dt$ is the material (Lagrangian) derivative with respect to time, and other variables are as defined previously. Without phase change or chemical mass transport $D\rho /Dt=0$. This results in:

$$\frac{D}{Dt}\left(\rho {\mathit{u}}_a\right)=\rho \frac{D{\mathit{u}}_a}{Dt}.$$

(26)

The temporal change in kinetic energy results from pre-multiplying (26) with ${\mathit{u}}_a$:

$$\rho {\mathit{u}}_a·\frac{D{\mathit{u}}_a}{Dt}=\rho \frac{1}{2}\frac{D{\left|{\mathit{u}}_a\right|}^2}{Dt},$$

(27)

where $\left|{\mathit{u}}_a\right|=\sqrt{{\mathit{u}}_a^T{\mathit{u}}_a}$. Returning to (26), the right-hand side is expanded

$$\begin{array}{cc}{\displaystyle \rho \frac{D{\mathit{u}}_a}{Dt}}\hfill & {\displaystyle =\rho \frac{D(\mathit{u}+{\mathit{u}}_d)}{Dt}=\rho \frac{D\mathit{u}}{Dt}+\rho {\mathit{a}}_d}\hfill \\ & {\displaystyle =\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·(\rho \mathit{u}\otimes \mathit{u})+\rho {\mathit{a}}_d.}\hfill \end{array}$$

(28)

Equation (25) is now cast in the non-inertial reference frame by assuming that the contribution of $\lambda \nabla ·\mathit{u}$ to the gas is negligible as per [38]:

$$\left[\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·(\rho \mathit{u}\otimes \mathit{u})+\rho {\mathit{a}}_d\right]=-\nabla p+\nabla ·(\mu +{\mu }_t)\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)+\rho \mathit{g}+\sigma \kappa \nabla \alpha .$$

(29)

In this work, we compute the mechanical energy in a manner that is consistent with the discretization of the governing equations. This is achieved by pre-multiplying the discretized governing equations at each computational cell i with ${\mathit{u}}_{a_i}$ followed by integration over time:

$$\begin{array}{cc}\hfill & {\displaystyle \sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\left(\frac{D\rho \mathit{u}}{Dt}+\rho {\mathit{a}}_d\right)\mathrm{d}V\mathrm{d}t=}\hfill \\ {\displaystyle }\hfill & {\displaystyle \sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\oint }_{\partial {\mathsf{\Omega }}_i}\left(\mathbf{\tau }·{\mathit{n}}_{{\mathsf{\Omega }}_i}\right)\mathrm{d}\mathsf{\Omega }\mathrm{d}t+\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\rho \mathit{g}\mathrm{d}V\mathrm{d}t+\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\sigma \kappa \nabla \alpha \mathrm{d}V\mathrm{d}t,}\hfill \end{array}$$

(30)

Here, ${\mathit{u}}_{a_i}=\frac{1}{V_i}{\int }_{V_i}{\mathit{u}}_{a}dV$ is the volume averaged absolute velocity and ${\sum }_i{V}_i=V$. Further, ${\mathsf{\Omega }}_i$ and ${\mathit{n}}_{{\mathsf{\Omega }}_i}$ are, respectively, the surface of $V_i$ and outward pointing unit normal.

We can proceed to group the terms in (30) into various energy containers. The potential $E_P$ and surface $E_{\sigma }$ energies are, respectively, given by

$$E_P\left(t\right)=\sum _i{E}_{P,i}=-\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\rho \mathit{g}\mathrm{d}V\mathrm{d}t\mathrm{and}$$

(31)

$$E_{\sigma }\left(t\right)=\sum _i{E}_{\sigma ,i}=-\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\sigma \kappa \nabla \alpha \mathrm{d}V\mathrm{d}t.$$

(32)

The stress tensor is divided up into internal and external contributions. Internal contributions are inside the tank boundary and include compressive energy $E_C$ (work performed on the fluid through pressure), as well as viscous dissipation $E_V$. We write these for a finite volume $V_i$ (shown in Figure 3), and to obtain the total for the domain, we simply integrate over all cells in the domain:

$$E_C\left(t\right)=\sum _i{E}_{C,i}=\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\oint }_{\partial {\mathsf{\Omega }}_i}p{\mathit{n}}_{{\mathsf{\Omega }}_i}\mathrm{d}{\mathsf{\Omega }}_i\mathrm{d}t$$

(33)

$$\begin{array}{cc}{\displaystyle E_V\left(t\right)}\hfill & {\displaystyle =\sum _i{E}_{V,i}={\int }_0^t\left({\mathcal{P}}_V+{\mathcal{P}}_V^{turb}\right)\mathrm{d}t}\hfill \\ & {\displaystyle =\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\oint }_{\partial {\mathsf{\Omega }}_i}\left[(\mu +{\mu }_t)\left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right)·{\mathit{n}}_{{\mathsf{\Omega }}_i}\right]\mathrm{d}{\mathsf{\Omega }}_i\mathrm{d}t,}\hfill \end{array}$$

(34)

with ${\mathsf{\Omega }}_i$ any internal surface bounding internal volume $V_i$. The external energy addition, or work performed on the fluid by the domain surface $\mathsf{\Omega }$ or work performed on the fluid by non-inertial forces, $W_{ext}$, is given by

$$\begin{array}{ccc}{\displaystyle W_{ext}\left(t\right)}\hfill & =\hfill & {\displaystyle -{\int }_0^t{\mathit{u}}_d·{\oint }_{\partial \mathsf{\Omega }}p{\mathit{n}}_{\mathsf{\Omega }}\mathrm{d}\mathsf{\Omega }\mathrm{d}t+}\hfill \\ & \hfill & {\displaystyle +{\int }_0^t{\mathit{u}}_d·{\oint }_{\partial \mathsf{\Omega }}\left[(\mu +{\mu }_t)\left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right)·{\mathit{n}}_{\mathsf{\Omega }}\right]\mathrm{d}\mathsf{\Omega }\mathrm{d}t.}\hfill \end{array}$$

(35)

where, in this work, the domain boundary velocity ${\mathit{u}}_d$ is constant in space.

For a flow field that varies smoothly in time, from (27), the kinetic energy contained in $V_i$ at any time t is:

$$\begin{array}{cc}{\displaystyle {\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\left[\frac{\partial \rho \mathit{u}}{\partial t}+{\nabla }^T(\rho \mathit{u}\otimes \mathit{u})+\rho {\mathit{a}}_d\right]\mathrm{d}V\mathrm{d}t}\hfill & {\displaystyle ={\int }_0^t{\int }_{V_i}\rho \frac{1}{2}\frac{D{\left|{\mathit{u}}_a\right|}^2}{Dt}\mathrm{d}V\mathrm{d}t}\hfill \\ {\displaystyle }\hfill & {\displaystyle ={\int }_{V_i}\rho \frac{1}{2}{\left|{\mathit{u}}_a\right|}^2\mathrm{d}V+E_{K,0},}\hfill \end{array}$$

(36)

with $E_{K,0}$ the kinetic energy at $t=0$. However, liquid slosh involves impact events (liquid–liquid and liquid–solid), which results in an immediate loss in kinetic energy [39]. Marrone et al. [27] observed that this can be accounted for via the artificial compressibility scheme [40,41]. For the fractional step method employed in the present work, this irreversible energy loss ($E_{loss}$) is calculated as the difference between the left- and right-hand sides of (36) as per [22]:

$$E_{loss}\left(t\right)=E_K\left(t\right)-\sum _i{\int }_0^t{\mathit{u}}_{a_i}·{\int }_{V_i}\left[\frac{\partial \rho \mathit{u}}{\partial t}+{\nabla }^T(\rho \mathit{u}\otimes \mathit{u})+\rho {\mathit{a}}_d\right]\mathrm{d}V\mathrm{d}t.$$

(37)

where

$$E_K=\sum _i{\int }_{V_i}\rho \frac{1}{2}{\left|{\mathit{u}}_a\right|}^2\mathrm{d}V-E_{K,0}$$

(38)

The total irreversible dissipated energy is then

$$E_{diss}\left(t\right)=E_{loss}\left(t\right)+E_V\left(t\right)$$

(39)

We can write (30) at any node i and re-shuffle to cast in the following form for the mechanical energy computed at any time, assuming the initial energy is zero

$$W_{ext}=E_K+E_P+E_C+E_{\sigma }-E_{diss},$$

(40)

The terms denote (from left to right) the work performed on the fluid through the non-inertial boundary surface, kinetic energy, potential energy, compressive (or dilatation) work, surface energy and dissipation. Equation (40) states that the total energy added to the fluid through the domain boundary must equal the sum of the stated energy containers on the right-hand side of the equation or the total energy inside the fluid. In Elemental^®, the contribution of each energy container to the total energy balance in (40) is calculated consistently with the respective contribution from each term in the momentum conservation Equation (29). This is completed very carefully when the contributions from each term are added to the fractional step solution method. Since momentum is conserved, this ensures that (40) is satisfied and that all the energy added to the fluid through the accelerating boundary is accounted for.

5. Results

We now present simulation results from both SPH (liquid only) and finite volume methods (liquid and gas). For 2D simulations, both numerical methods used small perturbations of the initial conditions to perform a statistical analysis of some result metrics. A subsection will be dedicated for each method to report on two- and three-dimensional results for both liquids modeled viz. oil and water.

5.1. SPH-Flow Solver Results

SPH-flow results, analyzed in this section, are obtained in a two-dimensional framework and compared with 3D results, respectively, for both water and oil. As discussed in [12], the problem studied in this work exhibits a chaotic nature. To highlight this behaviour, ten different 2D simulations have been performed introducing in the initial rest condition a small random noise on the particles positions.

5.1.1. Two-Dimensional SPH Water Results

We start to consider the test with water. Figure 4 shows the free-surface particles just prior to the first roof impact for the first four runs. Clearly, at this stage, the interface position from four different solutions are starting to deviate. The spatial resolution adopted is $N:=H/\Delta x=$ 200, corresponding to about 133,000 particles in the liquid domain.

Figure 5 depicts the vorticity field of four simulations, doubling the spatial resolution from $N=$ 50 up to 400. The time instant considered is related to the fourth bottom liquid impact event. The chaotic nature of the flow is also clearly visible in this plot, where the free-surface of the four fluid domains present large differences from each other. Furthermore, it is quite evident in the case of water that, by increasing the resolution, the number of eddies increases considerably. This is an indication that we are still resolving the inertial range of the turbulence regime. This is also confirmed by the fact that for the first three resolutions (N = 50, 100 and 200), the turbulent viscosity is one order of magnitude greater than the material or laminar viscosity. In a more in depth-analysis shown in [12,13], a resolution of $N=800$ is used for obtaining good LES simulations (i.e., for resolving the Taylor microscales).

However, because of this chaotic behavior, ten runs are performed for the four different spatial resolutions (N = 50, 100, 200 and 400) in order to perform a statistical analysis for controlling the convergence of the numerical results. On the other hand, the CPU costs for resolution N = 800 is too large to perform ten runs, and even for a 2D framework, a cluster machine would be required. The limitation is not mainly linked to the number of particles, which are about two million for N = 800, but to the fact that more than 1.5 million time iterations are required. Figure 6 reports the ensemble average of the dissipated energy of the forty simulations performed for the spatial resolution N = 50, 100, 200, 400. A convergence of the curves is obtained, where the error bars refer to the standard deviation measured between the set of the ten simulations. Even if a resolution N = 800 would be required, the result trend for N = 200 and 400 indicates that the average dissipation will likely remain within the standard deviation of the N = 400 case. It is worth noting that the results for the oil test case presented here are slightly different with respect to the ones reported in [13] since a different acceleration profile was used for the tank motion.

5.1.2. Two-Dimensional SPH Oil Results

Considering the oil case, the average dissipation results are much closer (compared to water), as shown in the plot of Figure 7. For all the three resolutions $N=50$, 100, 200, the averaged curves are quite close to each other, with the difference being less than the standard deviation measured. At the Reynolds number Re = 4660, even with the coarser resolution, it is possible to capture the main vorticity structure of the wall boundary layer as well as the large eddies inside the fluid. Even though the vortical structures in Figure 8 seem different at N = 400 compared to N = 50, the turbulent viscosity magnitude always remains of the same order as its laminar counterpart for all the simulations. Since the curves at N = 50, 100 and 200 are very close, the ensemble average is not performed at N = 400, for which only one simulation was performed.

5.1.3. Mechanisms of Fluid Dissipation

The plots in Figure 9 show the energy dissipation due to real and turbulent viscosity at three different spatial resolutions for both oil and water cases. Looking at the water (left of Figure 9), it is clear that the viscous dissipation due to the sub-scale model is predominant compared to the real viscosity. This is expected, considering the high Reynolds number. However, it is interesting that the absolute value of $E_V^{turb}$ does not increase with N while for the real viscosity component, it does. This indicates that the viscous dissipation related to the real viscosity can be resolved directly beyond the sub-scale turbulent model, which is something expected when a good LES is performed. For the oil cases, $E_V^{turb}$ largely decreases with N and remains lower than $E_V$, indicating that the LES simulation is much closer to being resolved, compared to water. This is again expected with the large difference in Reynolds number between the two liquid cases.

Figure 10 depicts the energy dissipated by the numerical diffusion of the scheme, which acts mainly during the liquid impact events. Because of the complexities of the flow studied, this component—like the $E_V^{turb}$ component—is not negligible, even for the lower Reynolds oil case.

In Figure 11, the time history of the energy dissipation $E_{diss}$ for the water test-cases and for one of the runs at resolution $N=400$ is compared with the histories of the mechanical energy (kinetic plus potential energy) in the Ni-FoR, along with $W_{NF}$ and $W_{dyn}^{ext}$. As already commented in Section 2.2, the mechanical energy of the fluid in the Ni-FoR remains quite limited and small during the experiments with respect to the dissipated energy. Furthermore, in this plot, it can be seen how both the work containers $W_{NF}$ and $W_{dyn}^{ext}$ tend to $E_{diss}$ at the end of the experiment $t_f=30$ T, indeed already after fifteen periods they are practically the same.

Figure 12 depicts the comparison between the time histories of the tank acceleration $a_{tank}$, the vertical velocity of the fluid center of mass ${\dot{y}}_G\left(t\right)$ and the power exerted by the non-inertial forces on the fluid $P_{NF}$. The latter is characterized by large peaks when the first two are in counter phase giving an indication of when intense fluid dissipation mechanisms are taking place.

5.1.4. Three-Dimensional SPH Water Results

Figure 13 shows the 3D free surface for the water test-case at the first roof impact. The 3D flow behavior is similar to the 2D one in the sense that during the first upward acceleration of the tank, a Rayleigh–Taylor instability is influenced by free-surface waves, which are generated due to the collapse of the menisci at the lateral walls. The waves move from the lateral walls towards the center of the tank. The further flow evolution is characterized by a complex 3D free surface, as well as a complex interaction of three-dimensional vortex tubes, which makes physics much more complex compared to the two-dimensional case. In addition to these differences, the dimensionless energy dissipation $E_{diss}/\Delta E$ remains quite close, as reported in

Table 1. Total energy fluid dissipation, $E_{diss}\left(t_f\right)$ for 2D/3D SPH simulations with water and with oil using different spatial resolution $N=H/\Delta x$. Numerical outputs are compared with the experimental data of [4]. The 2D results are ensemble averages of ten different simulations, while 3D simulations are given by single simulations.

Single-Phase $\mathit{\delta }$-LES-SPH Results (SPH-Flow Solver)
	Test-Case	N = 50	N = 100	N = 200	N = 400
$E_{diss}\left(t_f\right)/\Delta E_w$	2D test with water	−10.3	−12.4	−13.2	−13.6
$E_{diss}\left(t_f\right)/\Delta E_o$	2D test with oil	−14.0	−14.4	−14.6	−
		N = 22	N = 33	N = 50	N = 75	Exp. Data
$E_{diss}\left(t_f\right)/\Delta E_w$	3D test with water	−13.1	−14.5	−15.3	−15.9	−16.3
$E_{diss}\left(t_f\right)/\Delta E_o$	3D test with oil	−17.8	−17.7	−17.6	−17.9	−16.7

(for more details, see also [42]). We need to underline again that the 2D results are ensemble averages of ten simulations, while 3D is given by single simulations. Because of the CPU costs, the spatial resolutions adopted for the 3D simulations are $N=22$, 33, 50 and 75. As expected for the water-case, $E_{diss}$ varies significantly with N and, conversely to the 2D case, there is not a clear convergence trend in the results. In the water test-case, the highest resolution is still too coarse to obtain a convergent result and an effective LES simulation. A spatial resolution as high as $N=400$ corresponds to about 430 million particles and would require about one year of computing time on a cluster with 100 nodes. Even with the above constraints, the dissipated energy given by the $3D$ simulations at resolution $N=75$ is closest to the experimental data by [4].

5.1.5. Three-Dimensional SPH Oil Results

In the oil 3D test-case $E_{diss}$ is evaluated with four resolutions and—as for the 2D simulations—remains fairly close to each other. From Table 1, the oil 3D results overestimate the experimental data by approximately 7%, while 2D simulations undervalue the dissipation by approximately 13% and this is linked mainly to the missing friction effects on the lateral walls (for more details, see also [42]).

It is important to point out that the experimental uncertainty and its effect on $E_{diss}$ are unknown. The experimental value of $E_{diss}$ is calculated using Equation (12), where just small shifting errors between the forces and the tank velocity measures may imply non-negligible errors on the dissipation. For the above reasons, an in-depth investigation on the experimental side is still ongoing [11].

Figure 14 depicts the time histories of different energy components, including real, turbulent and numerical dissipation obtained for the maximum spatial resolution N = 75. Analogous to the 2D results, in the water case, the turbulent energy dissipation is the largest, while in the oil case, the dissipation due to the real viscosity is dominant.

In Section 3, the role of center of mass motion of the fluid is highlighted in the dissipated energy. The top plot of Figure 15 depicts the time histories of the vertical position $y_G\left(t\right)$ of the fluid center of mass for the water test-case evaluated by the SPH-Flow solver for both 2D and 3D framework. Remarkably, the two time histories are very close, up to $t=10$ T, while large differences appear in the rest of the evolution. However, in the first ten periods, about 85% of slosh dissipation takes place. In the bottom plot of Figure 15, the $y_G\left(t\right)$ time histories are reported for the oil test case. As already commented above, the differences in the flow behavior between 2D and 3D are larger here. In particular, the plots of Figure 15 show that the fluid center of mass in 3D is subjected to a larger motion, which is related to the higher slosh dissipation reported in Table 1.

5.2. Finite Volume Method Results with Elemental^®

5.2.1. Simulation Setup

For 2D and 3D simulations in Elemental^®, a structured, Cartesian mesh is used. The mesh resolutions used are similar to the particle resolution in the SPH simulations, as shown in

Table 2. Mesh resolutions for different simulations.

	Mesh Resolution:
Mesh Description and Dimension	${\mathit{N}}_{\mathit{x}}$×${\mathit{N}}_{\mathit{y}}$ (2D) or	Cell Dimension $\left(\mathsf{\mu }\mathbf{m}\right)$
	${\mathit{N}}_{\mathit{x}}$×${\mathit{N}}_{\mathit{y}}$× ${\mathit{N}}_{\mathit{z}}$ Nodes (3D)
$N50$ 2D	160 × 96	625.00
$N100$ 2D	320 × 192	312.50
$N200$ 2D	640 × 384	156.25
$N33$ 3D	125 × 75 × 75	800.00
$N50$ 3D	160 × 96 × 96	625.00

.

The interface is initialized as a flat, horizontal plane. To model the effect of liquid meniscii at the walls, liquid wedges are initialized as equilateral triangles with a height of $1.5$ mm, as explained in [12] and used in the SPH simulations. As mentioned earlier, interface perturbations are applied to the interface in 2D simulations. These perturbations are introduced as sinusoidal waves on the interface, all with an amplitude of $0.15$ mm, or $1/10$ the height of the meniscus. Five different wavelengths with randomized phase angles for each are used, resulting in differences in the interface shape once the Rayleigh–Taylor instability is formed with the first upward acceleration of the tank. Since the amplitudes of these perturbations were so small, the total difference in mass between simulations was negligible. Figure 16 shows the interface position for four different simulations at the $N50$ resolution for oil and water to illustrate the effect of the initial interface perturbations.

All simulations are run using the weakly compressible gas formulation for air and treating the liquid as an incompressible liquid. A continuum surface force (CSF) surface tension model is used to model surface tension. Like the SPH method, a Smagorinsky–Lilly large eddy simulation (LES) turbulence model is applied, with a Smagorinsky constant of 0.18 used throughout. For all simulations, the turbulent viscosity is capped to twice the laminar viscosity. The resulting linear system of Equation (23) is solved using a bi-conjugate gradient (BiCG) solver with an incomplete lower-upper (ILUT) preconditioner.

The measured tank acceleration is applied as a body force to all fluids in the tank, along with gravitational acceleration, exactly the same as in the single-phase approach. The evolution of the energy, as described in Section 4.4, is calculated by integrating (in time) the different energy components each time step. The total velocity at each time step is calculated as per Equation (14), where $\mathit{u}$ is the latest fluid velocity relative to the domain, and ${\mathit{u}}_d$ is the domain velocity as measured experimentally. Note that due to experimental measurement errors, there are slight inconsistencies between the measured velocity and acceleration of the tank. This can be seen, for example, in the “drift" of the potential energy. However, for the total dissipation, it can be shown that the inconsistency does not noticeably affect the calculated dissipation.

5.2.2. Oil 2D Results

In Figure 17, the oil interface is shown at different time instants as well as the vorticity in the fluid for a $N200$ resolution case. The first frame on the left at $t=1.20T$ shows the interface just after the central oil structures have impacted the roof for the first time. The asymmetry in the interface structure is a result of the randomized interface initialization explained earlier. Compared to the single-phase results in Figure 8, the liquid looks much more dispersed, or there is a higher degree of mixing between phases.

Since the spatial distribution in the tank differs from very early between single and two-phase results, it is difficult to perform a precise comparison of the vorticity. However, the overall vorticity magnitude in the liquid is similar between the approaches. Note that for the Elemental^® finite volume results, we use a scaled vorticity $\tilde{\mathbf{\omega }}$ throughout

$$\tilde{\mathbf{\omega }}=\frac{\mathbf{\omega }}{2\pi f_0}\frac{\mu }{{\mu }_{\ell }}=\frac{\nabla \times \mathit{u}}{2\pi f_0}\frac{\mu }{{\mu }_{\ell }}.$$

(41)

The last term is a viscosity ratio scaling factor, which is equal to 1 in the liquid phase and ${\mu }_g/{\mu }_{\ell }$ in the air. Without this scaling, the vorticity in the air phase is dominant and does not allow for a comparison with the SPH results in Figure 8.

To illustrate the effect of mesh resolution, the turbulent viscosity ratio is shown in Figure 18 at the same time for different mesh resolutions. The areas where the turbulent viscosity ratios are largest are mostly in the air. It was observed for all resolutions that there is large vorticity generated in the air phase as it is accelerated around in the tank by the liquid. As the resolution increases, it is evident that the inertial range for turbulence is increasingly resolved.

The energy distribution for a $N200$ case is represented in Figure 19, along with a magnification of the first five cycles. The total energy dissipation is shown in the light peach shading. In the zoomed plot (bottom), it can be seen that the dissipation starts to increase gradually after $t=1.20T$, which corresponds with the first roof impact of the oil. The sharp increase in dissipation just after $t=2.05T$ corresponds to the second roof impact. All the contributions to the total dissipated energy are presented in the plot, with the total viscous dissipation having the largest contribution. The total work done on the fluid, $W_{NF}$, consists of the work done by the pressure, $W_{NF}^P$ at the boundary, as well as the viscous work, $W_{NF}^V$. It is clear that the viscous component is negligible, and the majority of the work is performed by the pressure. The compressive $E_C$ and surface work $E_{\sigma }$ are not plotted here as both were negligible.

The impact-related loss term is significant, yet slightly less than the total viscous dissipation. The influence of mesh resolution on the total dissipation is shown in Figure 20. The difference between mesh resolutions is significantly more than in the SPH case, Figure 7, but the average dissipation has a convergent trend with increasing resolution. The individual dissipation components, as described by Equation (39), are presented in Figure 21. The impact loss is very similar for the $N50$ and $N200$ resolutions, yet the $N100$ resolution shows slightly more dissipation. There is no clear explanation for this difference and further investigation is required. Looking at the viscous dissipation, however, there is a clear trend with mesh resolution. The laminar component of the dissipation is not converged, while the turbulent component remains very similar between the three resolutions, indicating that most of the inertial range is resolved, as also shown earlier in Figure 18.

5.2.3. Water 2D Results

Figure 22 shows the location of the liquid phase at various stages, as well as the vorticity in the liquid. The same vorticity scaling is applied as for the oil. It was noted that for oil, the phases seemed more dispersed in comparison to the SPH results. For water, this tendency is similar. Comparing the 2D approach results of oil and water, the water contains even more small structures, such as droplets and ligaments. The vorticity magnitude is comparable to oil and the SPH results, but the vortical structures seem finer than in the oil, as was also observed in the SPH results.

When considering the turbulent viscosity ratio for water at different resolution levels, Figure 23, it is clear that this higher Reynolds number case has much more unresolved turbulence compared to oil, as expected.

Figure 24 shows the energy distribution for a 2D water simulation at the $N200$ resolution. Similar to the oil case, there is very little dissipation in the first cycle of the tank up to the first roof impact, whereafter there is a gradual increase with the first ($t\approx 1.20T$) and sharp increase with the second roof impact ($t\approx 2.05T$). Also similar is the negligible contribution of $W_{NF}^V$, $E_C$ and $E_{\sigma }$. What is different in the case of water is that the impact-loss-related dissipation is much more significant than the viscous dissipation.

The total dissipation is plotted in Figure 25, which again shows a larger variation between resolutions compared to the SPH results, but apparently convergent with increasing resolution.

The dissipation is divided into the impact loss and viscous components, plotted at the top and bottom of Figure 26. The impact-related losses are much more dominant, compared to oil and are also convergent. The viscous dissipation is much lower compared to the oil case. As can also be seen in Figure 23, the viscous dissipation is not resolved and both the laminar and turbulent components increase significantly with resolution. This can be expected since the Reynolds number is almost two orders of magnitude larger than the oil case.

5.2.4. Oil 3D Results

Three-dimensional Elemental^® runs were completed at resolutions of $N33$ and $N50$, as per Table 2. Figure 27 shows the oil–air interface in light brown, represented by an iso-surface of the VOF at the 0.5 level. The developing Rayleigh–Taylor instability is shown at $t=1.11T$, just before the first roof impact. The interface is very three-dimensional in nature, which explains the difference in computed energy dissipation as compared to the 2D simulations (see below). Different from the 2D simulations, no initial perturbations were imposed on the interface resulting in a symmetrical interface topology. After the first roof impact at $t\approx 1.17T$, the flow becomes very chaotic with large mixing of the two phases.

Figure 28 shows the energy distribution for the 3D oil simulation at the $N50$ resolution, with similar information as contained in the 2D results in Figure 19. Similar to the 2D simulations, the viscous component of boundary work $W_{NF}^V$ is negligible. For the dissipation in this simulation, the impact loss component is much more significant than the viscous component, which is different from the much finer $N200$ resolution in 2D. A comparison of impact loss and viscous dissipation is presented in Figure 29, including results from 2D. Comparing between 2D and 3D simulations at the $N50$ resolution level, both the laminar and turbulent components are larger in 3D than in 2D, with the turbulent (LES) component significantly so. However, a finer resolution will be required when considering the $N200$ resolution in 2D, which still has a much more significant laminar viscous component. Impact losses are also larger in 3D than in 2D, as shown in the middle plot. The result is that the total dissipation at the $N50$ resolution in 3D is more than the finest resolution ($N200$) in 2D. This indicates that the problem is truly three-dimensional in nature.

5.2.5. Water 3D Results

An iso-surface of the water interface for an $N50$ resolution 3D simulation is shown in Figure 30. As with oil, the interface shape is very symmetrical before impact but also strongly three-dimensional. The interface perturbation from the static meniscus grows into sheets of liquid with larger horn-like structures in the center. These horns and small streams in the corners are first to impact the roof, just before $t=1.11T$. To illustrate how mixed the two phases become after the first rounds of roof and floor impact, three views of the interface at $t=2.21T$ are shown in the last row.

The energy distribution at the $N50$ resolution is shown in Figure 31. In this plot, it can be seen that the dissipation through impact losses is a dominant component of the total dissipation for this resolution. The impact loss is clearly a much larger component in the total dissipation when compared to the 3D oil result at the same resolution, as shown in Figure 28. For viscous dissipation, the turbulent component is dominant.

Figure 32 presents a comparison between 3D resolutions, as well as 2D for water dissipation components. Similar to the oil 2D vs. 3D comparisons (for the two-phase model), the 3D simulations show clearly higher dissipation than in 2D. For impact losses, even the very coarse 3D simulation ($N33$) has higher impact-related dissipation than the $N200$ 2D simulation.

Table 3. Total energy fluid dissipation, $E_{diss}\left(t_f\right)$ for 2D/3D Elemental^® simulations at spatial resolution $N=H/\Delta x$. Numerical outputs are compared with the experimental data of [4]. The 2D results are ensemble averages of ten different simulations, while 3D results are given by single simulations.

	Elemental® Total Dissipation Results
	2D			3D		Exp. Data
	N = 50	N = 100	N = 200	N = 33	N = 50	-
$E_{diss}\left(t_f\right)/\Delta E_w$ Water	−6.6	−9.2	−10.6	−9.3	−12.2	−16.7
$E_{diss}\left(t_f\right)/\Delta E_o$ Oil	−7.5	−9.8	−11.4	−9.4	−12.5	−16.3

shows the dissipation results for Elemental^® at the end of the simulation compared to experimental energy dissipation estimates. The 3D predicted values correlate more closely to the experimental values as compared to the 2D simulations. This is as expected due to the strongly 3D nature of the Rayleigh–Taylor instability. Although the 3D predicted dissipated energy is not yet converged, the results seem to trend toward the measured values with mesh refinement.

6. Discussion on the Comparison between Elemental and SPH-Flow Results

In the previous sections, it was underlined that the problem studied in this work involves complex three-dimension flows involving complex liquid–gas interactions. Furthermore, for the water test-case, the high Reynolds number involved requires very high spatial resolution in order to perform effective large eddy simulations. In a 3D context, the CPU costs of those simulations are still prohibitive for both the adopted solvers, Elemental^® and SPH-flow. As a consequence, in the present work for the water test-case, convergent results are not achieved, and the effect of the spatial resolution N is not negligible in the results presented. However, from the point of view of the slosh dissipation, the single-phase model and 2D framework are able to supply a reasonable estimation of $E_{diss}$. This result is in agreement with previous literature on different kinds of problems (see, e.g., [24,26,43,44]).

Regarding the SPH-Flow solver when increasing the spatial resolution N, the 3D results for both water and oil test-case show a tendency to overestimate the dissipation compared to experimental measurement. This may be linked to the action of the gas phase, which is neglected in the adopted $\delta$-LES-SPH model. Conversely, in the case of the Elemental^® solver, it is evident that the gas–liquid interface is more complex when using a two-phase model. The liquid and gas phases are much more dispersed compared to the single-phase model. Consequently, very high spatial resolutions are needed to resolve the smallest scales (i.e., the small gas bubbles, the small liquid drops and the thin liquid jets).

The comparison between the two solvers and, therefore, the effect of the gas–liquid interaction is confirmed in 2D results of the vertical position of the fluid center of mass $y_G\left(t\right)$, depicted in Figure 33. In these plots, the SPH-Flow solver manifests a significantly larger motion of the fluid center of mass for both water and oil test-cases. In the water test-case, as expected, the solvers present remarkable differences not only in the fluid amplitude but also in the signal time phases. This is not the case in the oil results, where the two solvers show a good agreement in the time phase of $y_G\left(t\right)$ and a reduced difference in the amplitude of liquid motion.

Figure 34 depicts the vertical motion of the fluid center of mass for the three-dimensional simulations of the two solvers. The result trends are similar to the 2D simulations. In fact, the center of mass motion of the SPH model presents a larger amplitude again, and in the water test-case, a significant time-shifting appears on the $y_G\left(t\right)$ curves in the comparison between the two codes. The oil case also presents a significant amplitude difference in center of mass motion, even though the motion is in phase with the SPH results.

CPU Costs

Regarding the CPU costs of the two solvers for the test-cases with water:

• The SPH-Flow 3D simulation with N = 50 (750,000 particles) requires 36 h on 192 cores for 225,000 time iterations ($\Delta t\simeq 18\mathsf{\mu }\mathrm{s}$).
• The 3D Elemental simulation with N = 50 (1,500,000 cells for discretizing both liquid and gas phases) requires 46 h on 192 cores for 225,000 time iterations (the same as SPH-Flow).

The SPH simulations ran on the “Liger” supercomputer at Ecole Centrale de Nantes, which is equipped with 12–core Intel Xeon (Haswell) E5–2680v3 processors. We evaluated a computational speed defined as:

$$\eta =\frac{CP{U}_{time}{N}_{cores}}{N_{iteration}{N}_{particles}}$$

(42)

For the SPH-Flow solver $\eta =148\mathsf{\mu }\mathrm{s}$, while for Elemental^® $\eta =100\mathsf{\mu }\mathrm{s}$; therefore, the performance of the two solvers are of the same order.

7. Conclusions

The violent sloshing flow studied in this work presents several complexities. While the SPH results are promising in terms of comparison against experimental data, the two-phase results suggest that there are gas–liquid interactions in the flow leading to enhanced mixing, which is not captured in the single-phase approach. Conversely, the physical model adopted in the Elemental^® solver is more rich, showing a more complex evolution of the gas–liquid interface compared to the single-phase model. However, the spatial resolution required for converged energy metrics seems to be significantly higher compared to the single-phase SPH model. This is clear when examining the dissipation components related to the viscous dissipation, including the sub-grid model. With the SPH solver, these dissipation components seem to remain bounded for the water test-cases and decreases for the oil test-cases with increasing resolution. Conversely, with the Elemental^® solver, the energy dissipation is not yet fully converged with increasing mesh resolution.

For a future and further step of what is presented in this work, it would be interesting to extend the analysis using a two-phase SPH model, as well as conduct simulations with Elemental^® at increased resolutions. Both these improvements, however, will require a large leap in computational resources.