LNG Tank Sloshing Simulation of Multidegree Motions Based on Modified 3D MPS Method

Modified 3D Moving Particle Semi-Implicit (MPS) method is used to complete the numerical simulation of the fluid sloshing in LNG tank under multidegree excitation motion, which is compared with the results of experiments and 2D calculations obtained by other scholars to verify the reliability. +e cubic spline kernel functions used in Smoothed Particle Hydrodynamics (SPH) method are adopted to reduce the deviation caused by consecutive two times weighted average calculations; the boundary conditions and the determination of free surface particles are modified to improve the computational stability and accuracy of 3D calculation. +e tank is under forced multidegree excitation motion to simulate the real conditions of LNG ships, the pressures and the free surfaces at different times are given to verify the accuracy of 3D simulation, and the free surface and the splashed particles can be simulated more exactly.


Introduction
For the LNG/LPG ships sailing in the waves on the sea, sloshing of the fluid in the tank is generated because of the six-degrees-of-freedom motions, which may induce destruction of inner membrane and corners. Sloshing is a complex nonlinear phenomenon; at first, the main approach is experiment, and then Abramson first proposed using linear potential theories to predict the fluid impact in the cylinder and sphere tanks [1]. Faltinsen and Timokha in the book "Sloshing" introduced the background and basis of the analytic method in detail to study the sloshing problems [2]. e analytic method is not appropriate to solve the nonlinear sloshing problem, so the CFD numerical method is applied to simulate the violent sloshing including the free surface whirling and breaking. e governing equation is Navier-Stokes equation, which is very time-consuming to be solved, and the effective simulation of the free surface changing over time is a main difficult problem.
ere are several ways of calculation method to simulate fluid motions and free surface; the first type is based on Eulerian grid system, the typical approaches are VOF (volume of Fluid) method [3] and level-set method [4], and fixed or adaptive grids are applied to describe the free surface.
e second type is based on the method of Euler-Lagrange system; the boundary and free surface are tracked by the Lagrange system, but the fluid field is solved by the Euler method; at present the widely used methods are Particle-In-Cell (PIC) [5] and Marker and Cell (MAC) [6]. e third type is meshfree method, which is based on Lagrange system; commonly used methods are Smoothed Particle Hydrodynamics (SPH) [7,8], Dissipative Particle dynamics (DPD) [9], MPS [10], meshless finite element [11], and Consistent Particle Method (CPM) [12]. Meshfree method is a kind of absolute Lagrange particle solution; since the divergence is avoided, the large deformation of the free surface can be simulated perfectly.
SPH was first used in astrophysics, and then it has been extended to be applied in hydromechanics and solid mechanics, which is developed by many scholars. Grenier et al. [13] presented a new SPH method to solve multi-fluid-flows problem based on Colagrossi and Landrini's theory [14], and the modification is related to the formulation proposed by Hu and Adams [15]. Agertz et al. [16] carried out comparison study on grid and SPH of interacting multiphase fluids; the differences are shown in detail through numerical computation examples. DJ Price discussed the treatment of discontinuities to solve Kelvin-Helmholtz instabilities across entropy gradients in SPH [17].
MPS is different from SPH in the numerical format, computational efficiency, and convergence; in SPH, the pressure is solved by state equations, and the time integral function is explicit; in MPS, the pressure is solved from implicit Poisson equation, and the time integral form is semi-implicit. In SPH, the kernel function is not only weight function, but the gradient is also used in calculation directly; in MPS, the kernel function is only weight function, and the gradient is completed with the aid of radial function. e basic principle of MPS is that a serial of particles are arranged in the fluid field to simulate the macroscopic fluid flow. e position, mass, momentum, and energy are obtained from the neighbor particles with the help of kernel function; in the Lagrange coordinates, the movement of every single particle is tracked to obtain the whole fluid field information. Koshizuka et al. applied MPS method to solve the breaking waves [10], vapor explosions [18], and jet breakup [19] problems.
Based on MPS, other models are invented; for instance, Gotoh and Fredso [20] proposed a multiphase flow model for solid-liquid simulation, Ikari et al. [21] proposed a gasliquid two-phase flow model, and Koshizuka model is modified to fit the phenomena treated in hydraulic engineering. Feng and Huang [22] presented a 3D MPS method to simulate the supersonic atomization process, and the stability enhancement is discussed in detail. e MPS model is widely applied to calculate the sloshing in tank problem, Zhang et al. [23] presented a comparative study of MPS and level-set method of sloshing flows simulation, Kim et al. [24] investigated the interaction between the vessel motion and the inner tank sloshing through CFD and MPS method, and Gong and Yao [25] used MPS to study the resonance sloshing phenomenon, and the advantage of MPS in simulating the free surface is demonstrated. However, nowadays, the researches are mainly focused on the 2dimensional simulation; although some flow problems can be solved by 2D method, 3-dimensional effect cannot be neglected; there is a certain degree of 3-dimensional effect on the free surface breaking and the fluid splashing, and in some practical engineering situations, the flow is a 3-dimensional problem, for example, the wave-making, the diffraction on the offshore structure, and the entry water problems; 3D sloshing simulation is widely developed in these years; A. Souto-Iglesias et al. [26] proposed a 3D SPH numerical method to solve antiroll tanks and sloshing type problems, Zhang et al. analyzed the second-order sloshing resonance phenomenon in a 3D tank based on the potential theory and perturbation techniques [27], and Liu et al. adopted a Youngs' VOF method based on finite difference method to simulate 3D large amplitude sloshing in a rectangular tank [28]; it is difficult to deal with these sloshing problems through 2D solution, so it is necessary to use 3D calculation to solve these flow problems. For the free surface particle detection, Tanaka and Masunaga [29] proposed a new representation of the incompressible condition to stabilize simulations and suppress the pressure oscillation; the free surface particles are detected more accurately in the pressure prediction process. Lee et al. applied the revised gradient model and improved free-surface-particle search method to simulate the free surface, and the numerical results are also compared against the experimental results [30]. A new geometric boundary detection technique for assemblies of spherical particles is described by Dilts [31]; the comparison with SPH on a ball-and-plate impact simulation shows qualitative improvement. Idelsohn et al. used Voronoi diagrams or Voronoi spheres to make the computation easier to recognize boundary surface nodes [32]. Oñate et al. discussed the application of particle finite element method for solving the problems in involving fluids with free surface and submerged or floating structures within a unified Lagrangian finite element framework [33].
Based on SPH solvers, Marrone S. et al. proposed a fast novel algorithm to detect the free surface in particle simulations, in which the accurate normal vectors to the free surface are made available [34]. Shibata et al. developed a new free surface boundary condition for the pressure Poisson equation by assuming that there are virtual particles over the surface [35]. Zheng et al. suggested a new free surface identification scheme, which is applied in the MLPG_R method; this method can improve the identification of free surface particles without increasing the CPU time [36]. Tsuruta et al. presented a new scheme with consideration of a potential in void space as space potential particle to reproduce physical motions of particles around free surface through a particle-void interaction [37]. In order to improve the accuracy and stability of the particle-based free surface flow simulation method, an alternative fluid interface particle detection technique with easy implementation, high accuracy, and low processing cost is proposed by Tsukamoto et al. [38]. Barecasco et al. used a sum of normalized relative position vectors from neighboring particles to test particle boundary status; the fluid domain is assumed to be covered by overlapping spheres [39].
In above references, the method is applied to predict the dam breaking model or sloshing free surface and pressures in rectangular tank; in this paper, the sloshing tank is a real LNG tank model, and free surface particle detection scheme is different from the previous method, and the multidegree excitation motions are considered in the prediction.

Numerical Calculation Method
where ρ is density of the fluid; t is time; u is velocity vector; p is pressure; ] is kinematic viscosity coefficient; F is body force.

Kernel Functions and Computational Operator Models.
In MPS method, the vector between two particles i and j can be obtained from the scalar quantities ϕ i and ϕ j at r i and r j , which can be expressed as a gradient form as follows: e Laplacian from is represented as where In the above equations, w(r) is kernel function, ϕ is arbitrary scalar function, d are the space dimensions, n 0 is the initial particle number density, and 〈n〉 i is particle number density at coordinate r i .
Kernel function is used to simulate the interaction between the particles; in MPS method, the differential operator is discretized by the weighted average of the particles, but the differential operator in SPH method is discretized on kernel function straightly, so the kernel function used in SPH must be second-order smooth.
In order to avoid accumulated deviation of two times weighted average calculation in the process of discretizing the divergence of shear stress tensor of momentum equation in traditional MPS method, the divergence discretization scheme of SPH method is implied in MPS to improve the computation accuracy; by this way, this modified MPS method is applicable for non-Newtonian fluid too.
In SPH method, the divergence discretization is where m is the mass of the particle, τ is shear stress tensor, and ρ i is the number density of particle i: For incompressible flow, the mass and the number density are constant, so the following divergence function can be obtained: where From the above equations, divergence operator is transformed into gradient calculation of kernel functions, so the accumulated deviation from weighted average calculation in every step is avoided; in equation (7), the cycle computing process brings about accumulated deviations, so when equation (8) is used instead of equation (7) the accumulated deviation can be avoided. Because equation (5) satisfies the conditions of SPH kernel functions, the cubic spline kernel functions used in SPH method can be adopted in MPS: where h is smoothing length.

3D Computation Extension of MPS Method. According to the pressure Poisson equation proposed by Khayyer and
Goth [40], the high-order spring terms should be considered; the equation can be written as In this numerical model, the above equation can be extended to 3D: Considering kernel equation (10), the above equation can be written as e 2D Laplace operator can be expended to 3D; for particle i and j, the corresponding scalar function is Mathematical Problems in Engineering Using the same treatment as the CPM method proposed by Luo et al. [41,42], according to Taylor expansion, ϕ jj can be written as erefore, neglecting the higher-order items, the following deduction can be derived: e three-dimensional form of higher-order Laplace operator is obtained eventually: According to the results of experiments, Abbas Khayyer and Hitoshi Goth [40] suggested adding a modified matrix to make the calculation more stable.
Based on Taylor expansion, p j can be decomposed at p i in 3D, which is presented by Feng and Huang as follows [22]: Multiplying ((w ij x ij /r 2 ij )dx) + ((w ij y ij /r 2 ij )dy) + ((w ij z ij /r 2 ij )dz) on each side of equation (18), Oger et al. proposed a normalization procedure to achieve the calculation [43], which can be used here to simplify the process; a 3D corrective matrix is introduced: 4 Mathematical Problems in Engineering where V ij is the statistical volume, which is defined as e normalized form can be written as

Artificial Viscosity Terms. A general equation proposed by Monaghan and
Kocharyan can be used in MPS method, which can be simplified as the following form [44]: In the above equation, F ij � (r ij /|r ij |) · ∇ i w ij , ρ ij is the average density between particles, and v sig is the maximum speed between particles; Price presented the following formulation to calculate this signal speed [17]: where C s is the sound speed and l 0 is the average distance between particles. Based on equations (19) and (25), equation (24) can be written as Given by 3D form,

Boundary Conditions.
In the numerical simulation, the boundary conditions should be modified in order to achieve the 3D calculation, which contain the free surface condition and boundary condition.

Neighbor Search Method.
In the Poisson equation, not all the particles are involved in the calculation; the pressure of the particles on the free surface should be assigned to zero so as to ensure that the results are not divergent. e traditional method is to identify the free surface particles according to the following discriminant: where n * i is the particle number density after modification; β is the coefficient of determination, which is less than 1.
e determination is according to the particle number density; this is not a very accurate method based on some researches of other scholars. So, in this paper, a neighbor particle searching method is proposed to determine the free surface; in MPS method, not only do the particles on the free surface should be identified, but also the particles close to the free surface should be dealt with as the free particles, so the searching distance is between 1 and 2.1 times of the particle distance; as shown in Figure 1, we assumed one particle located at the center of the sphere, the radius of the small sphere equals the distance between particles, and the radius of the big sphere is 2.1 times of the distance between particles. Based on the 2D neighbor searching method proposed by Pan et al. [45], which can be extended in 3D dimensional simulation, this space can be assumed to be divided into 48 parts, if there is no neighbor particle in 8 connected parts of this apace; then this particle can be considered to be on the free surface.
In order to improve computation efficiency, as Pan proposed [34], a number density criterion is given firstly: A low limit coefficient is added into equation (29); this is because when the number density is less than A, the particle is assumed to be at a broken region, which must be on the free surface obviously. rough this first step, the process can be simplified; not all of the particles are necessary to be judged as shown in Figure 1.
In the real computation, the neighbor searching process is very time-consuming, especially when the particles number is huge; in order to save the computing time, the above equation is inserted in the program to exclude the few and far particles. In addition, smashing and fragmentation may occur in the simulation process, and neighbors may appear in the adjacent areas of all particles, but it belongs to the case of free surface. At this time, the distribution of particles is sparse, and the density of particles is less than the general density of particles, which can be directly classified as free surface particles. rough these procedures, the neighbor searching is only operated for some of the particles; the computing time can be saved significantly.

Boundary Treatment.
In order to balance the internal pressure, two layers of particles are arranged close to the solid wall to simulate the boundary nonpressure condition, which are fixed at the definite location to avoid the fluid Mathematical Problems in Engineering particles flow out of the tank. Between the wall and the fluid particles, a boundary pressure layer of particles is arranged, which takes part in the pressure calculation, but the velocity and the displacement vectors are not changed, as shown in Figure 2(c).
If there is only one layer of boundary particles, the fluid particles close to the boundary may be determined as free surface particle, so other virtual nonpressure particles are arranged to the boundary; as shown in Figure 2, the blue particles are the boundary and the red ones are the virtual boundary pressure particles. In general, the density domain is 2.1l 0 , where l 0 is the initial distance between particles; if the range in equation (29) is too large, then another layer of boundary pressure particle is necessary to be arranged to avoid the wrong determination as free surface particle; this is proposed by Pan to simulate the 2D fluid sloshing in tank [46].

Collision of Particles.
e particles inside the fluid are colliding with each other in the kinematical motions; the repulsive force is introduced according to the local pressure, but it cannot be calculated when the collision happens around the free surface; the pressure of the particles on the free surface is set to zero, so the repulsive force cannot be generated when the particles are close. As Byung-Hyuk Lee (2011) mentioned [30], especially when the particles from outside collide with the free surface particles, the density is increased, which may not be recognized as free surface particle. In order to avoid this phenomenon, a collision model is employed; as shown in Figure 3, the initial distance between particles is l 0 , the velocity of the first particle is v 1 , the mass is m 1 , the velocity of the other particle is v 2 , and the mass is m 2 , when the distance is smaller than al 0 ; the repulsive velocity is defined as v ′ � − bv, where a and b are fixed coefficients.

Multidegree Excitation Motions
In general, the study of sloshing in tank is mainly to simulate the liquefied natural gas moving in the tank of LNG ship advancing in waves. In order to simulate the real situation of the tank sloshing in a LNG ship, two groups of multidegree motions are recommended to be calculated; from the longitudinal side, the pitching, surging, and heaving motions are coupled; from transverse side, the rolling, swaying, and heaving are coupled. e motions can be written as a unified equation: In the above equation, B i is the tank displacement at time t, A i is the amplitude of the i-th motion, T is the period, t is the time, and ϕ i is the initial phase.   Table 1.

Verification and Comparison of Modified Kernel
Function. In MPS method, there are many forms of kernel functions, which can be grossly divided into two types: one type is that when the distance is close to zero, the kernel function value is infinite, and, in the other type, when the distance is zero, the kernel function is equal to a finite value; the commonly used kernel functions are as follows: In (9), the smoothing lengthh � �������������� Δx 2 + Δy 2 + Δz 2 3 , where Δx, Δy, and Δz are the distances between particles; the kernel function curves are shown in Figure 7; KF1 is equation (31), KF2 is equation (32), and KF3 is equation (9). KF1 rises rapidly with the decrease of the distance between particles; KF2 and KF3 rise slowly. e value of KF1 is higher than KF2; KF3 crosses the curves of KF1 and KF2; when r is close to zero, the value of KF1 is unlimited. As shown in Figure 8, compared with the simulation results obtained by Yu using LEVEL-SET RANS method [48], which is similar to the experiment results, the simulation with KF1 is too dispersed, the liquid elevation of the simulation with KF2 is too high, and the simulation with KF3 is consistent with the level-set method.

Verification of Neighbor Particle Search Method.
In order to testify that the modified particle search method is effective in the numerical simulation, comparison between equations (28) and (29) is completed; when the search method is not effective, there would be particles in the internal fluid recognized as free surface particles; in this situation, the pressure must fluctuate very violently, so it can be seen from the predicted pressures which method is more effective. e calculations of the pressures of different cases by two methods are shown in Figure 9; using equation (25), there are few internal particles defined as free surface particles; this is the reason why the pressures are more smooth and reasonable.

Comparison and Simulation Results.
e free surfaces and the fluid with pressure distribution under the multidegree motions at different time steps are given in Figure 10; in order to show the sloshing clearly, the boundary and wall particles are not shown in the following 3D distribution figures.
e simulation results are compared with the simulations obtained by Yu [48], and the impact pressures  Figure 3: Collision model.  [47]. e pressures of the checking points under different conditions are shown in Figure 11; ten periods and one period predictions of every case are demonstrated in the following figures.
As shown in Figure 10, at the initial flow stage, the pressure is basically equal to static water pressure; when the water is slamming on the wall, the pressure of the reference points is increasing sharply, as shown in Figure 11, the pressures fluctuate with the time; there are two peak values of the pressure in some periods; the first is the impact of the fluid on the wall; after the first impact, some of the water goes up along the wall and then flows down under gravity, resulting in second impact on the bottom. e free surface is broken before impact on the wall; this is the reason why after the peak values the pressure fluctuates sharply. Ch1 and Ch11 are located on the inclined wall; as seen from Figure 11(a), the pressure goes up quickly and declines slowly, which agrees with the experiment values; the numerical predicted maximum values are underestimated compared with experiment data; a detailed one period figure is presented in Figure 11(a); although the predicted pressure curve is more smooth than the experiment curve, it is clear that the prediction is in high agreement with the experiment values.
Compared with case 1, the horizontal motion amplitude is larger and the vertical motion is smaller; the pitch is much smaller than that of case 1 too. Ch7 is located on the top of the tank; in this case, the 2D prediction cannot predict the slamming effect on the top; the 3D prediction of the pressure at Ch7, as seen from Figure 11(b), is consistent with experiment data. e difference between the 2D and 3D simulations is shown in Figure 10; it can be seen clearly that the sloshing on the top of the tank is obvious in 3D simulation. Because the upper tank is narrow, the particles are forced to move to the middle area; the pressure rises sharply when the free surface is higher than the knuckle of the inclined wall, which is the result of 3D calculation effect. Another obvious effect is that the free surface and the pressures are asymmetric between different periods, which is proven to be correct in the experiment. e 2D and 3D numerical predictions are presented in Figure 11(c); the point Ch11 is located below the lower knuckle of the upper chamfer; for case 3, the pressure of this point can be predicted by 2D method too; it is obviously seen Ch1 Ch7 Ch11 Figure 6: Tank simulation and pressure checking points.        [40]. A � amplitude of the excitation; A max � maximum amplitude of excitation � 50 mm; t � simulation time; ω � frequency of excitation � 2π/T; T �period of excitation � 1.5 s; d 0 � particle size � 3.0 mm. from Figure 10 that the free surface is slamming on the upper wall of the tank; at the same time, the pressure on the lower chamfer is obviously higher; the peak value of the pressure appears when the sloshing water is slamming on the checking point. e swirling motion in MPS simulation is presented as particles splashing phenomenon; the three-dimensional effect leads to the generation of the swirling flow.
Khayyer and Gotoh applied MPS-HS-HL-ECS-GC method to simulate a violent sloshing flow induced by periodic horizontal excitations [40]; the schematic sketch of the calculation is shown in Figure 12.
e simulation results are shown in Figure 13; as shown in above figures, the method proposed in this paper is distinguished by some unphysical irregularities in simulation as well as dispersed particles close to the free surface.
Based on the above results and analysis, compared with the 2D results, the pressure curves are similar; the reason is that the dimensional effect, which has great influence on the broken free surface and fluid splashing, is less important for the overall movement. ere is high frequency oscillation phenomenon in 2D calculation, but in 3D calculation the pressure curves are smoother; this is because the interaction scope is sphere; there are more neighbor particles, so the calculated pressure is smoother than the 2D results. e calculated 3D pressures are higher than the 2D results and closer to the experiment values. e pressure obtained in the experiment is negative in some time, which is caused by the air or air-water effect in the space between the sensor and the tank; for MPS simulation, the pressures are all positive.

Conclusions
e study shows that the 3D MPS method proposed in this paper can simulate the free surface and the broken surface; the splashed particles can be tracked. e initial positions of the particles are easily arranged; particularly for complex surface boundary problems, it is more convenient. Compared with 2D simulations and the results obtained by other scholars, the calculated pressures and the free surface are close and the 3D results are in high agreement with the experiments data; the 3D simulation can describe the detailed features of the sloshing, for instance, the broken free surface and the splashed particles. e comparison shows that 2D and 3D methods can both describe the overall fluid movement, so for the tank under single motion the 2D method may be more appropriate than the 3D method because of the higher computation quantity in 3D calculation, but when the tank is under multidegree motions and is not in a normal shape, the free surface and the splashed particles can be described more exactly by the 3D method, and for the problems of three-dimensional flow or other high requirements for the smoothness of the pressure field, 3D model simulation is more appropriate and effective.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.