Experimental Validation of Single- and Two-Phase Smoothed Particle Hydrodynamics on Sloshing in a Prismatic Tank

Abstract

This study aimed to validate the single-phase and two-phase smoothed particle hydrodynamics (SPH) on sloshing in a tank. There have been many studies on sloshing in tanks based on meshless particle methods, but few researchers have used a large number of particles because there is a limitation on the total number of particles when using only CPUs. Additionally, few studies have investigated the influence of air phase on tank sloshing based on two-phase SPH. In this study, a dedicated sloshing experiment was conducted at the National Research Institute of Fishing Engineering using a prismatic tank with a four-degrees-of-freedom forced oscillation machine. Three pressure gauges were used to measure local pressure near the corners of the tank. The sloshing experiment was repeated for two different filling ratios, amplitudes, and frequencies of external oscillation. Next, a GPU-accelerated three-dimensional SPH simulation of sloshing was performed using the same conditions as the experiment with a large number of particles. Lastly, two-dimensional sloshing simulations based on single-phase and two-phase SPH were carried out to determine the importance of the air phase in terms of tank sloshing. Based on systematic comparisons of the single-phase SPH, two-phase SPH, and experimental results, this paper presents a detailed discussion of the role of air-phase in terms of sloshing. The currently achievable accuracy when using SPH is demonstrated together with a few sensitivity analyses of SPH parameters.

1. Introduction

The increasing demand for liquefied natural gas (LNG) has had a significant influence on the capacity of LNG carriers. During the marine transportation of LNG, there is a dangerous phenomenon called sloshing that can be defined as the resonance of fluid inside a tank caused by external oscillations (i.e., ship motion in the case of marine transportation). When a ship is subjected to external oscillation that is close to the natural frequency of the tanks it is carrying, severe sloshing can occur and cause damage to tanks based on the violent movement of fluid, which leads to a high impact pressure on tank walls. This high impact pressure can cause explosions when volatile fluids, such as oil, are transported on a ship. In recent decades, the demand for sloshing analysis has increased significantly in the shipping industry. Because sloshing is defined by nonlinear free surface flows, both physical experiments and numerical methods are commonly used for such an analysis.

Computational fluid dynamics (CFD) has been widely adopted as an advanced numerical approach. Mesh-based CFD is a standard method that numerically approximates a space using meshes. There have been several pioneering works on sloshing simulation using mesh-based CFD. Kishev et al. [1] developed a novel CFD simulation approach based on constraint interpolation profiles. Kim [2] investigated both 2D and 3D liquid sloshing in containers by using the synchronized overlap-and-add algorithm. Recently, the OpenFOAM software was used by Jiang et al. [3] to simulate sloshing coupled with ship motions. Sloshing with different filling levels was investigated experimentally by Chen and Xue [4].

Recently, particle methods, which are also known as mesh-free CFD, have been applied to tank sloshing problems. Because particle methods are fully Lagrangian and truly mesh-free approaches, mesh handling is not necessary and a significantly deformed free surface can be captured intuitively without any special treatment. Smoothed particle hydrodynamics (SPH) [5,6] is a popular particle method that was first extended for free surface flows by Monaghan [7]. Since then, there have been many applications of SPH to free surface flow problems. Dalrymple and Rogers [8] demonstrated that SPH is able to model breaking waves on a beach in both two and three dimensions. Makris et al. [9] used SPH to model wave breaking over a relatively mild slope in two dimensions. Antuono et al. [10] used delta-SPH for water waves and drew comparisons to a boundary element method mixed-Eulerian-Lagrangian solver. Altomare et al. [11] used SPH to model water waves and performed verification using the second-order Stokes theory. By using a general stabilizing algorithm based on Fick’s law, Lind et al. [12] simulated wave propagation accurately without decay. Colagrossi et al. [13] investigated the dissipation mechanisms of gravity waves and found that the number of neighboring particles is an important parameter for the convergence of results. Recently, an experimental validation of SPH for long-distance propagation was carried out by Trimulyono and Hashimoto using a large numerical wave tank [14]. The application of SPH to fluid-structure interactions was reported by Antoci et al. [15], who modeled the interactions between fluids and flexible gates based on 2D simulations. Crespo et al. [16], Altomare et al. [17], and Barreiro et al. [18] used SPH to predict impact forces on coastal structures. Wei et al. [19] as well as Sarfaraz and Park [20] performed SPH simulations to investigate the impact of tsunamis on bridge piers. Marrone et al. [21] presented SPH simulations for modelling viscous flows around blunt bodies at low and moderate Reynolds numbers. Zisis et al. [22] executed multiphase SPH simulations by applying a robust SPH number-density scheme and compared their method to the arbitrary Lagrangian Eulerian method. Kawamura et al. [23] reported the validation of SPH for ship motion in waves. Gómez-Gesteira et al. [24] performed SPH simulations to capture green water overtopping with a fixed horizontal deck. Le Touzé et al. [25] applied SPH to green water and ship flooding problems. Recently, a versatile algorithm for open boundaries in SPH was developed by Tafuni et al. [26] for 2D and 3D complex flows such as water waves and flows past a surface-piercing extraterrestrial submarine. Gonzalez-Cao et al. [27] made comparisons between SPH and the volume of fluid method to assess violent wave collisions on coastal structures. Verbrugghe et al. [28] used the open boundaries to perform water wave generation and absorption in one wavelength of wave flume and validation was carried out using a theoretical solution.

The application of SPH to sloshing problems has been carried out by Delorme et al. [29], Landrini et al. [30], and Chen et al. [31] to validate impact pressures based on SPH with forced roll motion in two dimensions. Kim [32] used SPH and the finite difference method to validate impact pressures using a rectangular tank. Colagrossi et al. [33] corrected pressure fields as a solution to pressure noise using an improved SPH method and Shao et al. [34] used kernel gradient correction and a novel treatment of solid boundaries to achieve a smoother distribution of the pressure field for a rectangular tank. Rafiee et al. [35] performed 2D and 3D SPH simulations of a rectangular tank under sway motion using an improved SPH method. Bouscasse et al. [36], and De Chowdhury and Sannasiraj [37] investigated shallow water sloshing using delta-SPH under sway motion with both small and large amplitudes of excitation. Recently, long-duration simulations of sloshing at small filling ratios were performed by Green and Peiro [38], who found that the SPH results matched experimental results in terms of both surface elevation and force on the tank. These results clearly demonstrate that SPH is suitable for evaluating the sloshing of nonlinear free surface flows. However, most past studies used simple rectangular tanks and single-phase SPH methods with relatively small numbers of particles. In this study, sloshing in a prismatic tank with sway and roll oscillations was investigated and experimental validation of both single-phase and two-phase SPH methods was attempted.

To perform large-scale SPH simulations, general-purpose computing on graphics processing units (GPGPU) technology can be advantageous because it facilitates the use of many compute unified device architecture (CUDA) cores on GPUs. In this study, a large-scale 3D sloshing simulation was carried out using an open-source SPH solver called DualSPHysics [39] that can be downloaded at https://dual.sphysics.org. DualSPHsysics is based on the SPHysics code and consists of a set of C++ and CUDA libraries for handling real-world engineering problems, such as sloshing and water waves. GPGPU has been implemented in DualSPHysics [40], but multiGPU computation has not yet been implemented. Computation decreases with the implementation of multiple GPUs, as shown by Dominquez et al. [41]. There are three main steps for SPH computation on a GPU: generation of a neighbor list, computation of the forces between particles, and updating of physical quantities prior to the next iteration. These steps are all executed on the GPU until the end of simulation. Some initial processing is performed on the CPU, then all data are transferred to GPU, where they continue to evolve in the GPU memory. Only a small amount of data are occasionally passed back to the host to be written to the disk [40]. All simulations in this study were performed using CUDA cores on a single GPU to handle a large number of particles and reduce computation time. In SPH simulations, the number of available particles is crucial not only for accuracy, but also for the application to realistic ocean engineering problems. GPGPU technology can realize efficient parallel computing for fully explicit time-marching schemes, such as SPH, which results in reliable and feasible solutions.

To validate the GPGPU SPH simulations of sloshing phenomena, a dedicated model experiment was conducted at the National Research Institute of Fishing Engineering (NRIFE). External sway and roll motions were generated for a prismatic tank by a four-degrees-of-freedom oscillation machine. The sloshing experiment was repeated for several different filling ratios, amplitudes, and frequencies of forced oscillation. The fluid behavior in the tank and the hydrodynamic pressures at different points were recorded. Next, 3D sloshing simulations based on single-phase SPH accelerated by a GPU with large numbers of particles were carried out and the results were compared to the experimental results. We found certain limits on single-phase SPH in terms of reproducible accuracy for the impact pressure. Next, 2D GPGPU simulations based on two-phase SPH were carried out to analyze the effects of the air phase on sloshing. Two-phase SPH simulations were carried out for two filling ratios with external sway and roll oscillations. By applying two-phase SPH to sloshing, the prediction accuracy for both global fluid behavior and local impact pressure can be improved. Through a systematic comparison study, validation of single-phase and two-phase SPH for sloshing phenomena is carried out and the importance of the air phase in sloshing problems is elucidated.

2. Numerical Methodology

SPH was initially developed for astrophysical problems and it has been widely used in marine engineering since its first application to a free-surface flow by Monaghan [7]. SPH is a fully Lagrangian meshless method that adopts an interpolation scheme to approximate the physical values and derivatives of a continuous field using discrete evaluation points. These evaluation points are identified as particles that contain mass, velocity, and position values. These quantities are obtained as weighted averages of adjacent particles within a smoothing length (h) to reduce the range of contributions from remote particles. The main features of the SPH method, which is based on integral interpolants, are described in detail in References [42,43]. In SPH, the field function A(r) in a domain Ω can be approximated using the integral approximation in Equation (1), where W is the kernel function and r is a position vector.

$$A\left(\mathit{r}\right)={{\displaystyle \int }}_{\Omega }A\left({\mathit{r}}^{\prime }\right)W\left(\mathit{r}-{\mathit{r}}^{\prime },h\right)d{\mathit{r}}^{\prime }$$

(1)

Equation (1) can be approximated in a discrete form by replacing the integral with a summation over neighboring particles in the compact support of particle a at the spatial position r. Therefore, the particle approximation of Equation (1) becomes the following.

$$A\left({\mathit{r}}_a\right)\approx {{\displaystyle \sum }}_{b}A\left({\mathit{r}}_b\right)W\left({\mathit{r}}_a-{\mathit{r}}_b,h\right)\frac{m_b}{{\rho }_b},$$

(2)

$$A_a={{\displaystyle \sum }}_b{m}_b\frac{A_b}{{\rho }_b}{W}_{ab},$$

(3)

where the volume associated with neighboring particle b is m_b/ρ_b, where m, ρ denote mass and density, respectively. In discrete form, Equation (2) leads to an approximation of the function at particle a, as shown in Equation (3), where $W_{ab}=W\left({\mathit{r}}_a-{\mathit{r}}_b,h\right)$ is a kernel function. Therefore, the derivative of these interpolants can be expressed by Equation (4).

$$\nabla A_a\left({\mathit{r}}_a\right)={\displaystyle \sum }_b{m}_b\frac{A_b}{{\rho }_b}\nabla W_{ab}$$

(4)

The Lagrangian forms of the governing equations of the Navier-Stokes equation are

$$\frac{D\rho }{Dt}=-\rho \nabla \mathit{v},$$

(5)

$$\frac{Dv}{Dt}=-\frac{1}{\rho }\nabla P+\mathrm{g}+\mathit{\tau },$$

(6)

$$\frac{Dr}{Dt}=\mathit{v},$$

(7)

where D/Dt is the material derivative, ρ is density, P is pressure, $\mathit{v}$ is velocity, $\mathrm{g}$ is gravitational acceleration, and $\mathit{\tau }$ represents diffusion terms. The capability of the SPH model depends on the choice of a kernel function. Kernel approximation is one of the main sources of errors in the SPH method. The selection of a kernel function not only affects computational efficiency, but also determines the accuracy and stability of the SPH method, as demonstrated by Cao et al. [44]. The kernel function must satisfy certain conditions, including the normalization condition, compact support, positivity, monotonic decrease with an increase in distance symmetric property, delta function property, and smoothness [43]. In this study, the Wendland kernel function [45] was used because it is computationally efficient and sufficiently smooth, as mentioned by Cao et al. [44]. This function is shown in Equation (8).

$$W\left(q\right)={\alpha }_D{\left(1-\frac{q}{2}\right)}^4\left(2q+1\right) 0\le q\le 2,$$

(8)

where ${\alpha }_D$ is equal to (7/4) πh² in 2D and (21/164) πh³ in 3D. q is the non-dimensional distance between particles a and b that is defined by r/h. The momentum equation to be solved in SPH is defined in Equation (9) for a water phase and Equation (10) for an air phase. The two-phase SPH method in DualSPHysics is based on a study by Colagrossi and Landrini [46], and it was developed by Mokos et al. [47].

$$\frac{d{\mathit{v}}_a}{dt}=-{{\displaystyle \sum }}_b{m}_b\left(\frac{P_{a+P_b}}{{\rho }_a\cdot {\rho }_b}+{\Pi }_{ab}\right){\nabla }_a{W}_{ab}+\mathrm{g},$$

(9)

$$\frac{\mathit{d}{\mathit{v}}_{\mathit{a}}}{dt}=-{{\displaystyle \sum }}_b{m}_b\left(\frac{P_{a+P_b}}{{\rho }_a\cdot {\rho }_b}+{\Pi }_{ab}\right){\nabla }_a{W}_{ab}-2a{\rho }_a^2{{\displaystyle \sum }}_b\frac{m_b}{{\rho }_b}{\nabla }_a{W}_{ab},$$

(10)

$$\mathrm{where}{\Pi }_{ab}=\{\begin{array}{ll}\frac{-\alpha \overline{c_{ab}}{\mu }_{ab}}{\overline{{\rho }_{ab}}}& {\mathit{v}}_{ab}\cdot {\mathit{r}}_{ab}<0\\ 0& {\mathit{v}}_{ab}\cdot {\mathit{r}}_{ab}>0\end{array},$$

where P_a and P_b are the pressures of particles a and b, respectively. ${\Pi }_{ab}$ is an artificial viscosity term defined by Monaghan [48], where ${\mu }_{ab}=\frac{h{v}_{ab}·r_a}{r_{ab}^2+0.01h^2}$. Additionally, ${\mathit{v}}_{ab}={\mathit{v}}_a-{\mathit{v}}_b$, where r and v are position and velocity vectors, respectively, and $\overline{c_{ab}}=0.5\left(c_a+c_b\right)$ is the mean speed of sound. Additionally, α is a coefficient that needs to be tuned to achieve proper dissipation.

The additional term in delta-SPH (${\delta }_{\mathsf{\Phi }}$), which was defined by Molteni and Colagrossi [49], is detailed in the paper by Antuono et al. [50,51]. This term was introduced to suppress density fluctuations. The equation of continuity with additional term of delta-SPH is defined in Equation (11).

$$\frac{d{\rho }_a}{dt}={\displaystyle \sum }_b{m}_b{\mathit{v}}_{ab}\cdot {\nabla }_a{W}_{ab}+2{\delta }_{\mathsf{\Phi }}h{c}_0{\displaystyle \sum }_b\left({\rho }_{b-}{\rho }_a\right)\frac{{\mathit{r}}_{ab}\cdot {\nabla }_a{W}_{ab}}{{\mathit{r}}_{ab}^2}\frac{m_b}{{\rho }_b}$$

(11)

DualSPHysics is based on weakly compressible smoothed particle hydrodynamics (WCSPH), where fluids are treated as weakly compressible. Therefore, the state equation defined in Equation (12) is used to determine the pressure field based on the computed particle density, where $b=\frac{c_0^2{\rho }_0}{\gamma }$, c₀ = c(ρ₀) = $\sqrt{\frac{\partial P}{\partial \rho }}{|}_{\rho 0 }$, and γ = 7. c₀, ρ₀, and γ are the speed of sound at the reference density, the reference density itself, and the polytropic constant, respectively. In the two-phase SPH method, a modified state equation is introduced, as shown in Equation (13), where $a=1.5\mathrm{g}\left(\frac{{\rho }_w}{{\rho }_a}\right)L$. ρ_w, ρ_a, L, and X are the initial water density, air density, characteristic length of the domain, and constant background pressure, respectively.

$$P=b\left[{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }-1\right]$$

(12)

$$P=b\left[{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }-1\right]+X-a{\rho }^2$$

(13)

The time-stepping scheme used in this study is the second-order explicit symplectic defined by Leimkuhler et al. [52]. This is an integration algorithm that is accurate with a time complexity of O($\Delta t$). It utilizes the predictor and corrector stages defined below.

$${\mathit{r}}_a^{n+\frac{1}{2}}={\mathit{r}}_a^n+\frac{\Delta t}{2}{\mathit{v}}_a^n$$

(14)

$${\rho }_a^{n+\frac{1}{2}}={\rho }_a^n+\frac{\Delta t}{2}{D}_a^n$$

(15)

$${{\mathit{v}}_a}^{n+1}={v_a}^{n+\frac{1}{2}}+\frac{\Delta t}{2}{\mathit{F}}_a^{n+1/2}$$

(16)

$${\mathit{r}}_a^{n+\frac{1}{2}}={\mathit{r}}_a^{n+\frac{1}{2}}+\frac{\Delta t}{2}{\mathit{v}}_a^{n+1}$$

(17)

F consists of a pressure gradient term, diffusion term, and gravitational term, where the n-step indicates the time step and $\Delta t$ is the time increment of the n-step. D_a is the divergence. The time step is calculated using the technique proposed by Monaghan et al. [53] as follows.

$$\Delta t_f=CFL\cdot min\left(\Delta t_f,\Delta t_{cv}\right),$$

(18)

$$\Delta t_f=\underset{a}{min}\left(\sqrt{\frac{h}{\left|{\mathit{f}}_{\mathit{a}}\right|}}\right),$$

(19)

$$\Delta t_{cv}=\underset{a}{\min }\frac{h}{C_s+\underset{b}{max}\left|\frac{h{\mathit{v}}_{ab}\cdot {\mathit{r}}_{ab}}{\left({\mathit{r}}_{ab}^2+{\eta }^2\right)}\right|},$$

(20)

where $\Delta t_f$ is based on the force per unit mas ($\left|{\mathit{f}}_{\mathit{a}}\right|$) and Δt_cv combines the Courant and viscous time step controls. CFL is a coefficient in the range of 0.1 ≤ CFL ≤ 0.3 [38]. Dynamic boundary particles (DBPs) are adopted based on the work by Crespo et al. [54]. DBPs are boundary particles that satisfy the same equations as fluid particles, but they do not move according to the forces exerted on them. Instead, they either remain fixed in position or move according to an imposed/assigned motion function (i.e., moving objects, such as gates, wavemakers, or floating objects). The particle shifting algorithm [12] was used to handle anisotropic particle spacing. When using the shifting algorithm, particles are moved to an area with fewer particles, which allows the domain to maintain a uniform particle distribution and eliminates voids to prevent noisy results.

3. Experimental Conditions

An experiment dedicated to the sloshing problem was recently conducted at the NRIFE using a prismatic tank. Figure 1a presents an overview of a sloshing test using a four-degree of freedom forced oscillation machine. Four video cameras were used in the sloshing test. Two were GoPro Hero cameras (GoPro Inc., San Mateo, CA, USA) that were placed near the tank to record the water behavior inside the tank and the other two cameras were used to record the external motion of the tank. Figure 1b presents the prismatic tank, which resembles a membrane-type LNG tank. The tank is scaled to 1:125 and the length is 0.36 m, the width is 0.30 m, and the height is 0.21 m. In the experiment, the tank was subjected to external sway, heave, roll, and yaw oscillations (see

Table 1. Experimental conditions.

External Motion	Frequency (Hz)	Amplitude (mm/deg)
Sway	1.08	0.82, 1.63, 3.26, 6.52
	1.34
Roll	1.04	1.09, 2.17, 4.33, 8.66
	1.34
Heave	1.04	2.85, 5.7
	1.34
Yaw	1.08	3.0
	1.34

). Several amplitudes and frequencies were used for the forced oscillation test. The impact pressure caused by sloshing is an important factor because it can lead to structural damage. To measure the hydrodynamic pressure on the tank wall at different locations, pressure gauges (SSK Co., Ltd., Tokyo, JAPAN) were positioned at three locations (P1, P2, and P3), as shown in Figure 1b. P1 is near the bottom of the sidewall, and P2 and P3 are on the upper chamber. All gauges were positioned near corners to capture the impact pressure. Sloshing tests were conducted at filling ratios of 25% and 50%. The water height was 5.25 cm for the filling ratio of 25% and P1 was near the calm-water free surface. The water height was 10.5 cm for the filling ratio of 50%. A schematic view of the prismatic tank is presented in Figure 2.

Table 2. Dimensions of the prismatic tank (units in meters).

Height (H)	0.21
Width (l)	0.30
Length (L)	0.38
Water height (d)	0.0525 (25%)
-	0.1050 (50%)

lists the dimensions of the tank.

Figure 3 presents time histories of external oscillations with the maximum displacement imposed on the tank in each mode during the experiment. For all conditions, a simple linear ramping motion was applied for the first 10 seconds to avoid huge inertial forces. By directly imposing the measured displacements of the tank during SPH simulation, the SPH method can be validated in terms of fluid behavior, including complex free-surface evolution and the time history of local pressure (e.g., impact at the corners of the tank). In this paper, only the sway and roll motions are presented because sloshing cannot be observed when the tank is subjected to heave and yaw motions. Figure 4 presents the measured pressure at P1 with a filling ratio of 25%. In this paper, pressure refers to hydrodynamic pressure, which can be obtained by subtracting the analytical hydrostatic pressure under calm conditions from the measurement data. The sloshing caused by the roll motion is more violent when compared to that caused by the sway motion because the center of rotation is higher than the top of the tank for the roll motion.

The natural period of the prismatic tank was estimated using Equations (21) and (22), which were derived by Faltinsen and Timokha [55], where ω_n is the natural frequency of the i-mode for a rectangular tank, d represents the water height, and l represents the length of the free surface in the direction of tank movement. For a prismatic tank with a chamfered bottom, a correction factor was introduced in Equation (22), where δ₁ and δ₂ are the horizontal and vertical dimensions of the chamfer, respectively.

$${\omega }_n=\sqrt{\frac{i\pi gtanh\left(\frac{i\pi d}{l}\right)}{l}}$$

(21)

$$\frac{{{\omega }^{\prime }}_n^2}{{\omega }_n^2}=1-\frac{{\delta }_1{\delta }_2^{-1}sin{h}^2\left(\frac{\pi i{\delta }_2}{l}\right)-{\delta }_1{\delta }_2^{-1}si{n}^2\left(\frac{\pi i{\delta }_1}{l}\right)}{\pi isinh\left(\frac{2\pi id}{l}\right)}$$

(22)

4. Numerical Setup

A tank with the same dimensions as that used in the experiment was used for SPH simulation to reproduce the sloshing phenomena observed in the experiment. A filling ratio of 25% was selected and the pressure gauge at P1 was located near the water surface. External oscillation was imposed based on a prescribed motion derived from the measured time histories of displacement during the experiment. In the simulation, an extra 5 seconds of simulation without external motion was included before oscillation began to settle the fluid particles.

Table 3. Computational settings for DualSPHysics.

Kernel Function	Wendland
Time step algorithm	Sympletic
Artificial viscosity with coefficient α	0.07
Coefsound (corresponding speed of sound (m/s))	60.0 (43.0)
Particle spacing (mm)	0.8
Coefh	1.2
Delta-SPH (δϕ)	0.1
Simulation time (s)	30.0

lists the numerical parameters/conditions used in the SPH simulation. The Wendland kernel function was selected for the simulation. In this scenario, “coefh” is defined as h/dp√3, where dp is the particle spacing, for three dimensions and h/dp√2 for two dimensions. The artificial viscosity coefficient α needs to be tuned to achieve proper dissipation. In this study, a value of 0.07 was used for all simulations based on trial-and-error testing. The parameter “coefsound” is defined as Cs/√(gd_swl), where Cs and d_swl are the speed of sound and the height of water under calm conditions, respectively. Because the speed of sound has a significant influence on pressure fields based on the state equation, a coefsound value of 60 was derived through numerical trials mentioned later. The CFL number is a coefficient for calculating the optimal time step to satisfy the Courant-Friedrichs-Lewy condition. A delta-SPH coefficient of 0.1 was used to reduce density fluctuation, as recommended by Molteni and Colagrossi [49].

Table 4. Particle spacing and total numbers of particles for the convergence study.

External Motion	Particle Spacing (mm)	Total Number of Particles
Sway	4.0	98,766
	3.2	179,271
	2.4	420,414
	1.6	1,376,755
	0.8	10,501,723

lists the numbers of particles used for a convergence study based on particle spacing. These particles were analyzed in a 3D simulation for a 25% filling ratio with external sway. The ratios of the width of the tank (l) to the particle spacing (dp), the inter-particle distance, were 75, 93.75, 125, 187.5, and 375 for the particle spacing values listed in Table 4. Using a finer particle spacing for simulation increases spatial resolution. It also reduces the gap between boundary particles and fluid particles, which enhances the prediction accuracy of impact pressure when using DBPs. When using DBPs in SPH simulation, a gap between fluid particles and boundary particles occurs, which is caused by an artificial force exerted to the boundary particles. It makes the point measurements by the pressure probe rather difficult to set on exact positions on the wall. Therefore, the gap between boundary and fluid particles needs to be considered in order to contain the typical pressure probe. The gap of 1.5h was used for this study as recommended by the DualSPHysics basic numerical setup. Figure 5 presents the results of the convergence study in terms of dynamic pressure under sway oscillations with A = 6.52 mm and f = 1.04 Hz. The SPH results using particle spacing values of 4.0, 3.2, and 2.4 mm cannot be used to evaluate the shape of pressure variation in each swing, while those using particle spacing values of 1.6 and 0.8 mm can. When comparing two results, the results for a spacing value of 1.6 mm exhibit a prominent spike of negative pressure that was not observed in the experiment. Additionally, the root mean square errors for particle spacing values of 1.6 mm and 0.8 mm are 119.33 and 56.83, respectively, which indicates that better agreement with the experimental results was achieved with a particle spacing value of 0.8 mm. Although even finer particle spacing could provide slightly better results, computation time would increase significantly. Therefore, a particle spacing of 0.8 mm will be assumed for the SPH simulations described hereafter.

5. Results and Discussion

5.1. 3D Single-Phase SPH Simulation

The speed of sound plays an important role in WCSPH. However, an artificially reduced value is commonly used to avoid excessive computational cost because an increase in speed requires a shorter time step, which increases total computation time. Using a reduced value for the speed of sound could result in less accurate results for impact pressure problems. Therefore, in this study, 3D single-phase SPH simulations using different values for the speed of sound were executed to determine the minimum adequate value. Figure 6a presents a pressure comparison for P1 under roll oscillations. The selected coefsound values were 20, 40, 60, and 80. Coefsound of 80 provides the best results, but increases computation time by four times compared to a coefsound of 20. To balance accuracy and computation time, a coefsound of 60 was selected in this study.

The natural frequencies of the prismatic tank with different filling ratios were calculated using Equations (21) and (22) (see

Table 5. Natural frequencies of the prismatic tank at different filling ratios.

Filling Ratio	Natural Frequency of the Prismatic Tank (Hz)
25%	1.10
50%	1.43

). To validate SPH for violent sloshing, the frequency of external oscillation should be close to the natural frequency of the tank. Therefore, values of f = 1.08 Hz and 1.04 Hz were selected for sway and roll oscillations, respectively, at a filling ratio of 25%. Additionally, the maximum amplitude was selected for the validation of each external motion. In this condition, the computation time for 30 seconds of simulation time was 347 hours.

Figure 7 presents a comparison of free surface evolutions under sway oscillations. In this moderate sloshing situation, single-phase SPH can accurately reproduce the resonance of the fluid and free surface deformation. Figure 8 presents a comparison of free surface evolutions under roll oscillations. With this type of oscillation, the fluid motion in the tank becomes more violent, as shown in Figure 8, and the water front touches the top of the tank. The SPH simulation accurately reproduces the run-up, breaking, and splashing phenomena during sloshing, but the water does not reach the top of the tank. This could be because the smoothing scheme and artificial viscosity used in SPH led to slightly greater dissipation compared to the real physics. Figure 9 and Figure 10 present comparisons of pressure variation over time at P1 between the SPH and experimental results. In these graphs, the arithmetic means of pressure and peak pressure are drawn together. The dashed line represents the arithmetic mean of pressure and the dotted line represents the mean of peak pressure under each oscillation type over nine cycles. The results reveal that the mean pressure and the mean of peak pressure are overestimated by SPH for both sway and roll oscillations, even though over 10 million particles were used. Additionally, one can see that the fluctuation in peak pressure in the SPH results is more prominent than that in the experimental results.

5.2. 2D Two-Phase SPH Simulation

As described in the previous subsection, a 3D single-phase SPH method was applied to the prediction of sloshing in a prismatic tank. As a result, the violent resonance of fluid with complex free surface excitation caused by external sway and roll was qualitatively reproduced with an accurate impact pressure near the tank’s corners. However, the simulated results did diverge from the experimental results, especially the pressure peak values, even while using over 10M particles. To identify possible elements for reducing these discrepancies, the two-phase SPH simulation was performed under the same conditions to analyze the influence of air phase on sloshing in a prismatic tank. The application of two-phase SPH was reported by Mokos et al. [47,56] and a ratio of the speeds of sound in liquid and gas phases of 7.5 had good agreement with experimental results [56]. In this study, the speed of sound of water was determined to maintain a ratio above 7.0 when considering the real speed of sound in air (i.e., 343 m/s). For the sake of a fair comparison and a clear discussion of the effects of the air phase, a 2D simulation using a single-phase SPH method with the same particle spacing as the previous simulation was performed. The computation time for two-phase SPH simulation is much greater than for the single-phase SPH simulation because air particles must be calculated in addition to water particles. Additionally, the minimum time step satisfying the Mach number becomes smaller based on the greater speed of sound. As a result, it is infeasible to run 3D two-phase SPH for sloshing simulations with a large number of particles. However, it is not suitable to use a different particle spacing value and parameters for two-phase SPH when comparing the results to those of single-phase SPH. Therefore, we decided to run 2D sloshing simulations based on two-phase SPH using a consistent particle spacing of 0.8 mm. This is a reasonable choice because pure sway and pure roll motions were introduced in the experiment, which means the flows in the tank can be regarded as 2D, rather than 3D. These assumptions led to a reduction in the total number of particles to 90,057 for 2D two-phase SPH simulation at a filling ratio of 25%.

Table 6. Computational settings in DualSPHysics.

Kernel Function	Wendland
Time step algorithm	Sympletic
Artificial viscosity coefficient α	0.07
Coefsound for water and air (corresponding speed of sound (m/s))	65 and 478 (46.5 and 343.0)
Particle spacing (mm)	0.8
Coefh	0.95
CFL number	0.2
Delta-SPH (δφ)	0.1
Simulation time (seconds)	30.0

lists the parameters used for two-phase SPH simulation. As mentioned previously, the speed of sound has a significant influence on impact pressure and the ratio of the speeds of sound between air and water is important from the perspective of numerical stability. Figure 11 presents comparisons of the free surface evolutions between single-phase and two-phase SPH, as well as the experimental results. The computation time for 30 s of simulation time was 81 hours for two-phase and 9 hours for single-phase SPH simulations, respectively. In the results for two-phase simulation, the red-colored region is the air phase and the blue-colored region is the water phase. Two-phase SPH exhibits a smoother free surface shape compared to single-phase SPH. This could be because the number of neighboring particles near the free surface is increased with the inclusion of air particles. The two-phase simulation shows better agreement with the experimental results, especially in terms of the run-up along the side wall. The single-phase simulation tends to exhibit less run-up and the fluid easily detaches from the side wall based on the exclusion of an air phase. Figure 12 presents the pressure contours under static and dynamic conditions for a filling ratio of 25% under roll oscillation. One can see that both SPH simulations yield smooth hydrostatic pressure contours [14], but the two-phase SPH exhibits a smoother pressure distribution under dynamic conditions. In SPH, flow is simulated by solving partial differential equations locally with a compact support. However, the pressure field can be naturally smoothed over the entire liquid domain. Figure 13 presents a comparison of local pressure at P1 for sway and roll oscillations. The peak values tend to be overestimated by both single-phase and two-phase SPH. The mean value is also overestimated for sway, but good agreement can be observed for the roll. In general, the mean peak pressure is slightly more accurate for two-phase SPH. It is known that pressure fields generally exhibit oscillation based on the effects of density fluctuation in WCSPH, but pressure fluctuation is also suppressed by taking the air phase into account. The shape of the pressure curve, especially the tail of the pressure curve after the peak, shows better agreement with the experimental results for roll for two-phase SPH. It should be noted that the single-phase 2D SPH results show better agreement with the experimental results compared to the single-phase 3D SPH results in Figure 9 and Figure 10, even though the same particle spacing is used. This implies that the accuracy of 3D simulation tends to be worse than that of 2D simulation. The SPH method adopts a Lagrangian formulation, so a uniform distribution over time is an important element for the accurate analysis of flows. This uniform distribution might be more difficult to maintain in 3D simulations when compared to 2D simulations.

Figure 14 presents images of free surface evolution under roll oscillation with f =1.30 Hz and A = 8.66° for a filling ratio of 50%. For two-phase SPH, the front shape of the run-up can be maintained along the side wall and touches the top wall, which is similar to the experimental observations. In contrast, for single-phase SPH, the fluid easily detaches from the side wall, similar to the filling ratio of 25%. Figure 15 presents the pressure contours under static and dynamic conditions. A smoother pressure distribution can again be observed in the liquid phase for two-phase SPH and one can see that the presence of air clearly affects the behavior of liquids, especially in terms of run-up and spraying. Figure 16 presents a comparison of the local pressures at P1 and P2 under sway oscillation with f = 1.34 Hz and A = 6.52 mm at a filling ratio of 50%. In the two-phase SPH results, the arithmetic mean of pressure shows reasonable agreement with the experimental results and the mean of peak pressure also shows reasonable agreement at P1, whereas single-phase SPH failed to predict these values accurately. The prediction accuracy of pressure at P2 is poor because the pressure at this point is relatively spikey when compared to that at P1. However, the two-phase SPH results show much better agreement with the experimental results than the single-phase SPH results at P2. A similar trend of improvement in terms of pressure accuracy can be observed for the roll oscillation shown in Figure 17. The accuracy of peak pressure at P2 is worse than that in the sway case because the flow inside the tank is more violent and the impact pressure is much greater than that in the sway case. However, the importance of the air phase is clearly demonstrated because single-phase SPH cannot capture the impact pressure at all.

Figure 18 presents comparisons between the two-phase SPH results and experimental results under roll oscillation. The intension for this comparison is to highlight the air and water mixing that occurs during violent sloshing. In the two-phase SPH simulation, the existence of air pockets and an air-water mixture is captured, which matches the experiment. This is one of the major advantages of two-phase flow simulations in terms of deep discussion of the sloshing phenomenon itself. Mokos et al. [56] detailed the issue of air mixture creation and void creation in a two-phase flow. It is clear that much finer particle spacing is required to capture the creation process of air-water mixtures quantitatively. Further studies for the detailed analysis of air effects using a large number of particles are expected. In addition, for the engineering standpoints, it is expected to couple two-phase SPH with a mathematical model for ship motion to carry out fully coupled ship-tank simulations, which is similar to those discussed in References [57,58,59,60].

6. Conclusions

In this paper, single and two-phase SPH simulations accelerated by a GPU were performed to analyze sloshing in a tank. For quantitative validation, a dedicated sloshing experiment was conducted using a realistically shaped prismatic tank that was subjected to external sway and roll oscillations. The experiment was successful, and fluid behavior and impact pressure at three different locations were measured. Next, single-phase 3D SPH simulations were performed with several particle distances. It was confirmed that 0.8 mm of particle spacing is necessary to capture the shape of pressure variation. It was also confirmed that the speed of sound has a significant influence on the impact pressure and pressure oscillation. By comparing the simulation results to the experimental results, it was demonstrated that large-scale single-phase SPH, using over 10 million particles, fails to represent free surface evolution and impact pressure accurately. Next, a 2D two-phase simulation was performed and validated by comparing its results to the experimental results, as well as those of the two-dimensional single-phase SPH simulation. Based on comparisons under different conditions, the importance of an air phase in sloshing phenomena was clarified and improvement of impact pressure simulation accuracy near the tank corners in terms of mean and peak pressures was confirmed.