Feng-ze Xie, Wei-wen Zhao, De-cheng Wan

Computational Marine Hydrodynamics Lab (CMHL), School of Naval Architecture, Ocean and Civil Engineering,State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

Abstract: For violent sloshing, the flow field becomes complicated and 3-D effect is non-negligible. In addition to the excitation direction, the wave can also propagate perpendicular to the excitation direction. Due to the superposition of waves from different directions, the impact pressure imposed on the wall of the tank may increase. In this paper, our in-house solver MPSGPU-SJTU based on moving particle semi-implicit (MPS) method coupled with GPU techniques is employed for the liquid sloshing simulation,to study the factors leading to the 3-D effect. Firstly, a series of sloshing simulations are carried out to validate the reliability of present solver. Then, the sensitivity of 3-D effect against some parameters, such as excitation frequency, dimensions of the tank and filling ratio, is checked through numerical simulations. Time histories of pressure obtained by 2-D and 3-D simulations are compared to judge the occurrence of 3-D effect. It concludes that effects of those parameters are all significant.

Key words: Liquid sloshing, 3-D effect, sensitive parameters, moving particle semi-implicit (MPS) method, MPSGPU-SJTU solver

Introduction

Due to the uneven distribution of energy in different regions, a large amount of energy needs to be transported from one area to another every year.Liquid sloshing is a significant issue in the transportation of liquefied natural gas, oil and liquefied petroleum gas. The liquid inside a partially filled tank will oscillate violently and exert large impact pressure on the tank under the external excitation, especially when the excitation frequency is close to the natural frequency of sloshing. Therefore,many researchers have investigated the characters and mechanisms of liquid sloshing.

Many approaches have been used to study the impact loads induced by the sloshing. Kim et al.[1]applied marker and cell (MAC) method for simulating sloshing in 2-D and 3-D prismatic tanks and studied the influence of internal baffles on pressure impact.Lee et al.[2]adopted Flow3D solver based on the finite volume method (FVM) to simulate the sloshing in a 2-D tank and investigated the effects of turbulence,viscosity, density ratio and compressibility of liquid on the sloshing loads. Panigrahy et al.[3]conducted a series of experiments to investigate the pressure distributions in tanks without or with baffles. Graczyk and Moan[4]experimentally studied the pressure impact in the LNG tank with different filling ratios.Belakroum et al.[5]used finite element method (FEM)to simulate the sloshing of rectangular tank and predicted damping effect of different baffles on sloshing in tanks partially filled with liquid. Chen and Wu[6]used a time-independent finite difference method (FDM) to study the sloshing in the rectangular tank under coupled motion and investigated the effect of excitation angle on the longitudinal force. Ali Goudarzi and Reza Sabbagh-Yazdi[7]employed the volume of fluid (VOF) method to capture the free surface and discussed the nonlinear effects on liquid sloshing in rectangular tank. Song et al.[8]studied flow kinematics and impact pressure in a partially filled rectangular tank through a series of experiments. Cho and Kim[9]used experiments and the linear potential theory respectively to investigate the effect of porous baffles on sloshing reduction in the rectangular tank.Seo et al.[10]investigated the interaction between ship motion and internal sloshing by using Rankine panel method. Zhuang and Wan[11]used the naoe-FOAMSJTU solver based on FVM to conduct the numerical simulation of FPSO motion coupled with LNG sloshing and investigated the impact pressure on the bulkhead. Stephen et al.[12]investigated the interaction between the ship motion induced by the waves and the sloshing in the tank based on the time domain potential theory. However, some limitations and deficiencies of those methods mentioned above are also exposed. Experiments are time consuming and limited by space. It is hard for traditional mesh methods to handle the free surface with large deformation, especially fragmentations and splashing.

Particle method has the superiority to handle the problem of violent free surface[13]and has been successfully applied for simulating liquid sloshing in previous work. Delorme et al.[14]applied smoothed particle hydrodynamics (SPH) method to simulate the sloshing in the tank with shallow water and investigated impact pressure on the walls. Koh et al.[15]adopted consistent particle method (CPM) to investigate the suppression of constrained floating baffle to sloshing in LNG tank. Chen et al.[16]used SPH to simulate sloshing in a rolling rectangular tank and analysed the pressure curve of a point near the free surface. Wen and Wan[17]developed a solver based on multiphase MPS method and applied it to simulation of sloshing in the tank with two-layer or three-layer fluids. Chen et al.[18]combined MPS with FEM to simulate 3-D sloshing in the elastic tank and investigated variation of impact pressures with the change of Young’s modulus of the tank walls.

In this study, our in-house solver MPSGPU-SJTU based on the modified MPS is employed to simulate the sloshing. The paper is organized as follows: Firstly, the modified MPS method and the GPU solver are presented,respectively. Secondly, through comparisons between numerical and experimental results, the accuracy of present solver is verified in the simulation of swaying rectangular tank. Finally, a series of numerical simulations are conducted to show the effect of excitation frequency, length of the tank and filling ratio on the formation of 3-D effect.

1. MPS method

MPS method was proposed by Koshizuka and Oka[19]for viscous incompressible fluid. In this section,basic theories and discretization process of MPS is presented in detail.

1.1 Governing equations

The governing equations contain continuity equation and Navier-Stokes (N-S) equation.

where V , ρ , t , p , ν and g refer to the velocity vector of the particle, the fluid density, the time, the pressure, the kinematic viscosity and the acceleration vector of gravity, respectively.

1.2 Kernel function

In MPS method, the interaction between particles is controlled by kernel function, which plays a role as weight function in the discretization process. In order to avoid non-physical pressure oscillation, the kernel function presented by Zhang et al.[20]is employed here,that is:

where r ander refer to the distance between two particles and the radius of the particle interaction.

1.3 Gradient model

The gradient model is the weighted average of the gradient vectors of a physical quantity φ between particle i and all its neighboring particles j . In order to meet the law of conservation of momentum, the model proposed by Tanaka and Masunaga[21]is adopted as follows

where d is the number of space dimension, n0is the initial particle density, r is the position vector relative to the coordinate origin.

1.4 Divergence model

Divergence model is similar to the gradient model and it is used to discrete velocity divergence in the pressure poisson equation (PPE) as

1.5 Laplacian model

The Laplacian model is the weighted average of the distribution of a quantity φ from particle i to neighboring particle j , which is needed in the solution of viscosity term of N-S equation and the space discretization of the PPE.

where λ is a parameter, which is used to compensate the errors resulting from the limited range of the kernel function.

1.6 Pressure poisson equation

In MPS method, the distribution of pressure is obtained through solving PPE. In this paper, mixed source term method is used to solve PPE, which is developed by Tanaka and Masunaga[21].

1.7 Free surface detection

An improved surface particle detection method developed by Zhang et al.[20]is adopted in this paper,which can distinguish the particles of free surface from the others efficiently. This approach is based on the asymmetry distribution of neighbouring particles.

1.8 Boundary condition

Diagram of boundary particles is presented in Fig.1. There are multilayer particles arranged at the solid boundary. One layer of wall particles is arranged near the fluid particles and their pressure is obtained by solving PPE. Two layers of ghost particles are configured because fluid particles lack neighboring particles on the side of the solid wall. The pressure of ghost particles is obtained by extrapolation. Both the wall particles and the ghost particles won’t update their velocity and displacement after they gain the pressure.

Fig. 1 (Color online) Diagram of boundary particles

2. GPU acceleration technique

In our previous work[22], GPU acceleration technique has been introduced to the MPS method and an in-house solver MPSGPU-SJTU has been developed.

One chip area of GPU is in possession of more arithmetic logic units (ALU) than CPU. So GPU can achieve high floating point operations per second(FLOPS) and deal with large amounts of data at the same time[23]. MPSGPU-SJTU solver is developed by C language on the CUDA platform. CUDA is a general parallel computing platform released by NVIDIA, which enables GPU to solve complex computing problems. The CUDA platform consists of a host part and a device part. CPU is the host part while the GPU is the device part.

The flow chart of the MPSGPU-SJTU solver is shown in Fig. 2. At the beginning of the MPS simulation, the host part allocates internal memory for itself and video memory for the device part. Then, the initial data is transferred from the host part to the device part and the MPS code runs on the GPU. Once the option of data saving is triggered, the data will be copied from the device part to the host part. At the end of the MPS simulation, the memory is released by the host part.

Fig. 2 (Color online) MPS flow chart of GPU implementation

3. Validation

The accuracy and the reliability of MPSGPU-SJTU are verified in this section. The simulationconditionsarethe sameasthoseinthe experiment of Kangand Lee[24]. Thesketch ofTank A is shown in Fig. 3. Tank A is 0.35 m long ( LA), 0.80 m wide (WA), 0.50 m high ( HA) and the water is 0.15 m deep ( DA). The probe P0 is arranged at the middle of the wall perpendicular to the excitation direction and 0.115 m away from the bottom. The initial particlespacing ( dp ) is0.005m( DA/ dp =30) andthe time step is 5×10-4s.The density of fluid is 1 000 kg/m3. Tank A sways harmonically under the external excitation.

where A is the amplitude of motion ( A =0.02 m),ω is the excitation frequency (ω =4.967 rad/s).

Comparison of the pressure impact among numerical and experimental results is shown in Fig. 4.2-D and 3-D numerical results are both in good agreement with experimental results, which indicates that the 3-D effect doesn’t occur in this case.Add ition al n umer ical res ults wit h DA /d p being 25 and50arealsoshownintheFigure.Resultsof simulations with different resolutions maintain almost consistence, showing the stability of present solver.

Fig. 3 The sketch of tank A

Fig. 4 (Color online) Comparison of pressure time histories among experimental and numerical results

Figure 5 shows the free surface profiles in tanks.The numerical simulations of flow field agree well with the experimental results, even for highly non-linear phenomenon such as breaking waves.

As the numerical simulation has been reasonably verified against sloshing experimental results, the numerical tools are used to investigate the sensitivity of 3-D effect against various parameters.

4. Case study: 3-D effect

The 3-D effect of sloshing was obviously observed in our previous work[25], based on which, the further work is done in this section. The sketch of tank

Fig. 5 (Color online) Snapshots of the experimental results (first line), 2-D numerical results (second line) and 3-D numerical results (third line) for sloshing flows

B is shown in Fig. 6. Tank B is 0.94 m long ( LB)0.47 m wide (WB), 0.47 m height ( HB). The initial particle spacing is 0.005 m and the time step is 5×10-4s. Tanks sway harmonically under the external excitation according to Eq. (12). The amplitude of motion is set to 0.015 m. Nature frequency (ω0) of the fluid depending on the filling depth is given by

Fig. 6 The sketch of tank B

where W is the width of the rectangular tank, h is the depth of the water in the tank.

Three probes are arranged on the wall of tank B to measure the time histories of pressure and their specific locations are listed in Table 1. Table 2 summarizes a selection of numerical cases considered in this sloshing study.

Table 1 Arrangements of pressure probes

Table 2 A selection of sloshing cases

4.1 Sensitivity against excitation frequency

In order to check the sensitivity against the excitation frequency, tank B is moved harmonically with different frequencies.

The snapshots of flow fields at different time instants are presented in Figs. 7 and 8. 3-D effect can be observed obviously in the tank B, which is under the external excitation of nature frequency. It can be noticed that the flow field in rectangular tank is asymmetrical. When the fluid flows to the right wall,the fluid in the middle of the wall flows faster than that on both sides, impacting the roof of the tank.When the fluid flows to the left wall, the fluid on both sides of the wall flows faster than that in the middle,reaching a higher position of the wall. The pressure varies tremendously along the longitudinal direction at the same level. There is no obvious 3-D effect observed when the excitation frequency is equal to 1.10ω .Besides, the distribution of flow field is symmetrical.Figure 9 shows the velocity vectors of free surface in some typical snapshots. In case 5, many velocity vectors have the component along the length of the tank. It can be noticed that the fluid transfers from the middle of the tank to both sides at that moment. In case 6, most vectors are parallel to the zx plane.

Fig. 7 (Color online) Snapshots of flow field in Tank B (Case 5: h / H B = 0.5, ω = ω0 , L =L B )

Fig. 8 (Color online) Snapshots of flow field in tank B (Case 6: h / H B = 0.5, ω= 1.1ω 0, L =L B )

Fig. 9 (Color online) Velocity vectors of free surface in tanks with different excitation frequencies

Figure 10 shows time histories of the impact pressure. When the excitation frequency is equal to natural frequency, it can be seen that there is great difference between 2-D and 3-D results. Before 13 s,3-D result is very close to the 2-D result, which implies that 3-D effect hasn’t been formed. From 13 s to 17 s, the 3-D effect is at initial stage and the peak pressure increases gradually. The double pressure peak in the 3-D tank disappears and is taken place by single pressure peak. The phase difference between 2-D and 3-D results becomes more apparent. This may be due to the fact that traveling waves are generated along the longitudinal direction and the waves from different directions are superimposed on each other.After 17s, the 3-D effect is at mature stage, the peak pressure measured in 3-D tank is much higher than that in 2-D tank, which may cause serious damage to the structure of the tank. 3-D effect can be also observed in Case 4. The formation of 3-D effect in Case 4 is earlier than that in case 5. Compared with case 5, phase of 3-D result is lag of that of 2-D result in case 4. For other cases, the general trends of the two signals are similar, which indicates that 3-D effect doesn’t occur in 3-D tank. Figure 11 shows the comparison of average pressure peaks between 2-D and 3-D results in different cases. These two signals both first increase then decrease with the increase of excitation frequency and the values of peaks occur near the nature frequency. When the excitation frequency is away from the nature frequency, the average pressure peaks measured in 3-D tank are close to that in 2-D tank. When the excitation frequency is equal to the nature frequency, the average pressure peak measured in 3-D tank is as much as about twice of that measured in 2-D tank.According to those numerical examples, the excitation frequency is a significant parameter on the 3-D effect. 3-D effect can be more easily observed when the excitation frequency is close to the natural frequency.

4.2 Sensitivity against tank length

In this subsection, the relationship between 3-D effect and tank length is discussed. The lengths of tank C and D are 0.470 m (0.5LB) and 0.235 m (0.25 LB)respectively. The other dimensions are consistent with tank B. Their excitation frequencies are all set to nature frequency. The probes in tank C and D are located at the same height as the probes in tank B and at the midpoint of the tank length.

Figures 12 and 13 show the flow fields at different moments. 3-D effect can also be observed in tank C. In the front of the tank C, the vertical motion range of the fluid near the right wall is larger than that near the left wall. In the back of the tank C, the vertical motion range of the fluid near the left wall is larger than that near the right wall. It seems that the tank moved diagonally instead of swaying. There is no obvious 3-D effect observed in tank D. Figure 14 shows the velocity vectors of free surface in some typical snapshots. In case 8, almost all vectors parallel with the plane on which the excitation force lies. In case 9, most of vectors near the right wall point towards the upper right corner in the front of the tank.In order to fill the void left by the fluid, the vectors near the left wall point towards the lower left corner in the front of the tank.

Fig. 10 (Color online) Comparison of pressure time histories between 2-D and 3-D results at P2 ( h / H B= 0.5, L =L B)

Fig. 11 (Color online) Comparison of maximum pressure between 2-D and 3-D results

The comparison of pressure time histories among tanks with different lengths is presented in Fig. 15.Time history of pressure measured in tank D has good agreement with the result calculated by 2D-MPS,which indicates that 3-D effect isn’t formed in the tank with the shortest length. The main reason is that the length of the tank restricts the wave generation in the longitudinal direction. Time histories of pressure measured in tank B and tank C are quite different from results of 2-D model in peak pressure and phase. 3-D effect observed in tank B is earlier than that observed in tank C. This implies that the formation of 3-D effect in the tank becomes earlier with the increase of the length.

Fig. 12 (Color online) Snapshots of flow field in tank C (Case 9: h / H B = 0.5, ω = ω0 , L= 0.5L B )

Fig. 13 (Color online) Snapshots of flow field in tank D (Case 8: h / H B = 0.5, ω = ω0 , L= 0.25L B )

It is noted that the length of the tank affects 3-D effect greatly. With the increase of the tank length, 3-D effect becomes more conspicuous in 3-D tank.

4.3 Sensitivity against filling ratio

The influence of filling ratio on 3-D effect is presented in this subsection. Tank B with different filling ratios is moved under nature frequencies and snapshots are shown in Figs. 16, 17. It can be noted that flow fields at different 2-D cross sections in the tank B with filling ratio of 0.25 are completely different, which confirms that 3-D effect has been formed. Even if 3-D effect can be observed in the tank B with filling ratio of 0.25 or 0.50, it is more difficult to find fluid motion laws when the filling ratio is low.It can be noticed that the free surface in case 10 is not as smooth as that in case 6, but 3-D effect isn’t formed in the tank B with high filling ratio. This is because that the free surface is closer to the roof in case 10 and the roof plays a role as a baffle, which prevents the further development of 3-D effect. Figure 18 shows the velocity vectors of free surface in some typical snapshots. In case 1, almost all vectors of the fluid,which haven’t reached the right wall, point to the middle of the right wall. In case 10, except for the vectors of splashing liquid, the y components of vectors of the free surface are relatively small.

Fig. 14 (Color online) Velocity vectors of free surface in tanks with different lengths

Fig. 15 (Color online) Comparison of pressure time histories at P2 among tanks with different lengths ( h / H B= 0.5,ω =ω0 )

Figure 19 shows the time histories of pressure in tank B with different filling ratios. When the filling ratio is low, the double peak changes into the single peak gradually after 17 s in 3-D tank while the double peak still exists in 2-D tank. The difference of phase between 2-D and 3-D results and the peak pressure of 3-D result increase over time. When the filling ratio is high, there is no big difference between the 2-D and the 3-D results.

It is found that 3-D effect is sensitive to filling ratios. 3-D effect can be easily captured when the filling ratio is low.

5. Conclusions

In this paper, the MPSGPU-SJTU solver based on modified MPS is employed to investigate the 3-D effect in the swaying rectangular tanks. The numerical results correlate well with experimental data, confirming the reliability of this solver. The influences of excitation frequency, length of the tank and filling ratio on the formation of 3-D effect are assessed.Based on the results of simulations, conclusions can be summarized as follows:(1) 3-D effect has great influences on the time histories of pressure measured in 3-D tanks. When the 3-D effect has been formed, there are great discrepancies in peak pressure and phase between 2-D and 3-D results.

(2) Pressure and velocity field in 3-D tank change dramatically when the 3-D effect occurs. The pressure varies tremendously along the longitudinal direction at the same level. The y components of velocity vectors become larger.

Fig. 16 (Color online) Snapshots of flow field in Tank B (Case 1: h / H B = 0.25, ω = ω0 , L =L B )

Fig. 17 (Color online) Snapshots of flow field in Tank B (Case 10: h / H B = 0.75, ω = ω0 , L =L B )

Fig. 18 (Coloronline)Velocityvectors of free surface in tanks withdifferentfillingratios

(3) The 3-D effect can be observed obviously when the excitation frequency is close to the natural frequency. The maximum pressure measured in 3-D tank is much greater than that in 2-D tank. The phase difference between 2-D and 3-D results is very large.When the excitation frequency stays away from the natural frequency, the pressure time histories of 3-D result are in good agreement with those of 2-D result.

(4) The length of the tank affects 3-D effect greatly. There is no 3-D effect in the shortest tank, but 3-D effect can be observed obviously in the medium tank and the longest tank. Besides, the formation of 3-D effect in the longest tank is earlier than that in the medium tank.

Fig. 19 (Color online) Comparison of pressure time histories between 2-D and 3-D results (Case 1: P1, Case 10: P3,ω =ω0 , L =L B )

(5) 3-D effect is sensitive to filling ratios. When the filling ratio is low, 3-D effect is easy to be formed.When the filling ratio is high, the roof plays a role as a baffle preventing the further development of 3-D effect.

Another in-house solver naoe-FOAM-SJTU[26]based on FVM has been widely applied in the ocean engineering field[27]. In the future, naoe-FOAM-SJTU will be coupled with MPSGPU-SJTU to simulate the motion of the ship with tanks in waves. The details of flow field around ship will be captured by The Eulerian method while the sloshing in the tanks will be evaluated by the Lagrangian method.