- Email: [email protected]

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Sloshing in a rectangular tank based on SPH simulation X.Y. Cao, F.R. Ming, A.M. Zhang ∗ College of Shipbuilding and Ocean Engineering, Harbin Engineering University, Harbin 150001, China

a r t i c l e

i n f o

Article history: Received 22 January 2014 Received in revised form 30 May 2014 Accepted 14 June 2014 Keywords: Sloshing Three dimensional SPH Bafﬂe Sloshing pressure

a b s t r a c t A ship would generate signiﬁcant sloshing when subjected to underwater explosion loads; the sloshing will reduce the ship’s stability and even cause the ship to capsize when coupled with internal liquid sloshing. It is of great signiﬁcance to research on the characteristics of the sloshing loads in a tank to improve the ship’s stability and security. The liquid sloshing in a tank is a complex process characterized with nonlinearity and strong randomness, and large amplitude sloshing is a great challenge for both theoretical models and numerical algorithms. Yet as a meshfree method, Smoothed particle hydrodynamics (SPH) has great advantages of solving such large deformation problems because of the nature of self-adaptiveness and Lagrangian. This paper focuses on the SPH core issues, such as the accuracy and the stability of the kernel function and boundary treatments. Firstly, the accuracy and computational stability of four common SPH kernel functions are simply investigated by two simple cases, and a more appropriate kernel function is selected. Secondly, the dummy particles and a novel boundary treatment considering the boundary motion are applied. Furthermore, the laws of impact pressure of the two-dimensional tank under forced rolling with different excitation frequencies and excitation angles are studied. Then, the inﬂuences of a bafﬂe for the liquid sloshing in a two-dimensional tank under forced surging are analyzed, and the action mechanisms of the bafﬂe are summarized. Finally, the coupled motion of swaying and surging for a three-dimensional tank is studied, which aims to lay a foundation for further study on the inﬂuence of sloshing loads on real ship motions. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction Underwater explosions may cause large-amplitude sloshing motions and even structural damage to ships. In such cases, the inﬂuent, together with the liquid inside the liquid tanks, would cause additional damages to ship structures due to the sloshing, which is a threat to the structural strength and safety of the ship. It will cause a great threat to the structural strength and ship safety, especially when the frequency of liquid sloshing approaches the natural frequency of the ship, which would affect the performance of shipboard equipment seriously and even cause the ship to capsize. Therefore, the studies of the sloshing characteristics and mechanisms of liquid tanks have great signiﬁcances for the ship stability and security. Liquid sloshing in a tank is a complex nonlinear motion involving free surface, splashing and complex coupling motion of liquid and the hull. Therefore, many scholars have focused on this classic problem.

∗ Corresponding author at: College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China. Tel.: +86 45182518443; fax: +86 451 82518296. E-mail addresses: [email protected], [email protected] (A.M. Zhang). 0141-1187/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.apor.2014.06.006

In theoretical researches, for example, the analyses based on linear potential ﬂow theory [1–5], relative accurate results can be obtained only for the sloshing tanks of simple shapes without internal structures. It is essential to mention that Faltinsen et al. [1–3] made great contributions to theoretical studies of liquid sloshing. Through experimental study, one can get a more accurate and reliable result. Many works about the model experiments of liquid sloshing in tanks have been done for different excitations, ﬁllings, compartments shapes and internal structures [6–9], and also liquid sloshing effects for ship motions [10]. However, most of the existing experimental studies are for multi-scale models rather than real ships. Based on the limitations of the theoretical and experimental researches, a variety of numerical methods are introduced to solve the problem of liquid sloshing. But most researches concentrated on the grid-based approach, for which special algorithms may be needed when tracking the motion of free surface. Taking the following studies for example, Liu et al. [11] simulated the liquid sloshing phenomenon in three-dimensional tanks with the ﬁnite difference method and the VOF method; Nagashima et al. [12] analyzed the characteristics of sloshing loads in a rotating tank with the ﬁnite element method combined with Level set method; Lee et al. [13] studied the inﬂuences of turbulence, viscosity, density and compressibility ratio of liquid on the impact pressure with the ﬁnite volume method, and the VOF method was introduced to track the free surface; Wu et al. [14] studied the inﬂuences

242

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

of a bafﬂe on the natural frequency of a two-dimensional liquid tank by the ﬁnite difference method and ﬁctitious cell method, etc. For the problems with free surface of large deformation, the meshfree method shows great advantages. The SPH method is a meshfree Lagrangian method, which can directly track free surface and material interface. This enables us to model the complex ﬂuid dynamics. There are two principal approaches in SPH method regarding ﬂuid compressibility, namely: incompressible smoothed particle hydrodynamics and weakly compressible smoothed particle hydrodynamics. Compared with the incompressible one, the weakly compressible procedure needs a relatively small time step, but it is easy to be paralleled and the boundary conditions of free surface are implicitly satisﬁed [15], and thus it appears to be more suitable for the free surface ﬂow problems. The obvious drawback of the weakly compressible procedure is unphysical noises in the pressure ﬁeld. To overcome the problem, many solutions have been proposed. Monaghan et al. [16] introduced an artiﬁcial viscosity term to reduce the numerical noise and stabilize the algorithm. Colagrossi et al. [17] corrected the density calculation using moving least squares (MLS) method, which is used in this work and presents some good results. Molteni and Colagrossi [18] and Antuono et al. [19] added diffusive terms in both the continuity and energy equations, and the beneﬁts and drawbacks can be found in [20]. Moreover, many modiﬁed SPH methods with Riemann solver [21–24] were introduced to reduce the numerical noises and avoid the use of artiﬁcial viscosity term, and further improve the accuracy and conserve the linear momentum, the total mass, and the total energy. However, the computational cost of Riemann solver is expensive compared with the weakly compressible SPH method especially for large particle numbers [21]. The SPH method has been applied in many ﬁelds, such as multiphase ﬂows [25], structures of large deformation [26], underwater explosion [27] and other aspects [28,29]. It has the advantage of solving the liquid sloshing problems and showing good results. As for two-dimensional tanks, Delorme et al. [30] analyzed the characteristics of impact pressure with a rectangular tank in rolling motion, Shao et al. [31] proposed a new coupling boundary method to simulate the rectangular tanks in the rolling or swaying motion, Marsh et al. [32] studied liquid sloshing in containers of different shapes, Colagrossi et al. [33] and Gotoh et al. [34,35] studied the violent sloshing phenomenon by using improved SPH method with MLS integral interpolators and enhanced incompressible SPH method, respectively; as for three-dimensional tanks, Vorobyev et al. [36] simulated the centralized sloshing in a ﬁxed cylinder tank and also the tank with obstacles inside, Raﬁee et al. [37] simulated sloshing in a rectangular tank whose width is far less than the length under a swaying excitation, Shen et al. [38] investigated sloshing motion in a barge in irregular waves, Bai et al. [39] studied the sloshing pressure and velocity in liquid natural gas (LNG) tankers and extended the tank with vertical bafﬂes in rolling motion. Through these studies, it can be seen that good results can be obtained with SPH method for two-dimensional and three-dimensional sloshing tanks, especially in a single degree of freedom. However, studies are still needed for the inﬂuence of inner tank structures on liquid sloshing and also the coupled motions of multi-degrees of freedom of the tank. The paper is organized as follows: ﬁrstly, the accuracy and stability of four commonly used kernel functions in SPH method are studied, and a new method for calculating the pressure of boundary particles is introduced; secondly, the inﬂuences of the rolling angles and the excitation frequencies on impact pressure are studied; thirdly, the inﬂuences of the bafﬂe on liquid sloshing in tanks in surging motion are investigated; ﬁnally, the coupled motions of swaying and surging for three-dimensional sloshing tanks are simulated. Through the comparison between the numerical and the

experimental results, the accuracy and reliability of the numerical algorithm are veriﬁed. 2. Theoretical background 2.1. Discretized SPH equations Due to the ﬁnite compressibility of real ﬂuid, the time step will be very small when the actual equation of state is applied. Therefore, the ﬂuid is usually regarded as weakly compressible in the actual calculation, and the pressure ﬁeld is obtained by solving the equation P = P(,e). In general, when the ﬂuid pressure is less than 1 GPa, the entropy effect on the pressure ﬁeld can be neglected [40]. The ﬂuid density is only a function of pressure. The SPH equations consist of the mass conservation equation, the energy conservation equation and the momentum conservation equation. When the pressure of the ﬂow ﬁeld is small, the ﬂuid is regarded as barotropic, and the energy has little effect on the pressure ﬁeld. Therefore, the energy equation is not solved. According to the kernel approximation and particle approximation in the SPH method, the density equation and the momentum equation [41] are expressed as follows: Di = −i Dt

(v j − v i ) · ∇ Wi j Vj

j

Dvi = g − mj Dt

Pi i 2

j

+

Pj j2

∇ Wi j + ˛h

j

(ci + cj )(v i − v j ) · r ij (i + j )(r 2ij + 0.01h2 )

)∇ Wi j

(1)

Pi = P(i )

where P, m, , c, v, r, g denote pressure, mass, density, speed of sound, velocity, coordinates and acceleration of gravity, respectively. rij = ri − rj . Besides, W is the kernel function, and the improved Gaussian kernel function is selected, see in Section 2.2. The subscripts i and j represent a pair of interacting particles. Coefﬁcient ˛ is set as 0.03 in the following simulations. From the above equations, one may draw that the density of particles is determined from the distribution of particles and the pressure ﬁeld. The selection of continuity equation is very important. The continuity equation used in this paper relates to the relative velocity of the neighboring particles in the support domain. Although the continuity equation cannot exactly guarantee the mass conservation, it can effectively reduce the calculation errors [42]. When dealing with hydrodynamic problems, usually an artiﬁcial viscous term is introduced into the momentum equation to ensure the stability of calculation [17], and the Monaghan form [43] is chosen. When particles approach each other, the performance of the artiﬁcial viscous term is a repelling force and the opposite on the contrary. The equation of state is obtained from Tait equation [16] by adiabatic approximation relations [44]: Pi = B

i

0

−1

(2)

where B = 0 c02 / and the speed of sound is linked with the ﬂuid medium through B; = 7, 0 = 1000 kg/m3 ; c0 represents the speed of sound which is set as more than 10 times of the maximum speed of the ﬂow ﬁeld. 2.2. Selection of kernel function Kernel approximation is one of the main sources of errors in the SPH method. The selection of kernel function not only affects the computation efﬁciency of SPH method, but also determines the accuracy and stability of the method. Especially, when the particle distribution is irregular, the calculation accuracy of the kernel function and its derivatives decay seriously, which has signiﬁcant

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

243

Table 1 Kernel functions [17,45,50]. Function name

Expression

Cubic spline kernel

W (r, h) = ˛d

4 − 6r 2 + 3r 3

0≤r<1 3 (2 − r) 1≤r<2 2≤r

0 4 (2 − r) (1 + 2r) 0 ≤ r < 2 W (r, h) = ˛d 2≤r 0

Wendland kernel (WC2)

5

5

for 2D and 3D

(3 − r) − 6(2 − r) + 15(1 − r) 5 5 (3 − r) − 6(2 − r) 5 (3 − r)

Quintic spline kernel

W (r, h) = ˛d

Improved Gaussian kernel

W (r, h) = ˛d (e−(q/h) − C) q ≤ ı

5

0≤r<1 1≤r<2 2≤r<3

2

In this paper the Wendland kernel 31 [50] is used and renormalized to obtain a compact support of 2h. The curves of the four kernel functions and their derivatives of the ﬁrst two orders are shown in Fig. 1. It is obvious that all the four kernel functions belong to the bell-shaped function, so they have similar properties in the support domain. According to the stress instability condition, it is easy for the kernel functions of this type to produce a stress instability phenomenon [51]. Regarding the size of support domain: the radii of the support domain of the cubic spline kernel and the Wendland kernel are only twice of the smoothing length. With the same smoothing length, the number of interacting particles for these two functions is less than that for the other two kernel functions, which is beneﬁcial to computational efﬁciency. Besides, the simulation results will be better if the functions are sufﬁciently smooth. As shown in the ﬁgure, the quintic spline kernel functions, the improved Gaussian kernel functions

effects on the results. Gaussian kernel function was introduced into the SPH method earlier, but the cubic spline kernel function [45] is more commonly used in the literatures; besides, the quintic spline kernel function [45] and the Wendland kernel functions [50] became popular in SPH simulations recently [46–48]. The speciﬁc forms of these four commonly used kernel functions are shown in Table 1. where, q = xi − xj and r = q/h; d is the number of dimension. ˛d is a variable coefﬁcient related to the smoothing length and spatial dimensions. The Gaussian kernel function is given in an improved form [17]. Compared with the classical Gaussian kernel, the improved Gaussian kernel is exactly zero at r = ı, and the difference between the two kernels is less than 0.04% [49]. As for the improved Gaussian kernel, ı is typically set equal to 3h, so C = e−9 , ˛ = [hd d/2 (1 − 10e−9 )]

−1

(3)

3

2 Cubic First derivative Seconnd derivative

1.5

0.5

0

0

-1

W/

1

W/

1

-0.5

-2

-1

-3

-1.5

-4

-2 -2

-1

0 r

1

Wendland First derivative Seconnd derivative

2

-5 -2

2

-1

(a) Cubic spline kernel

0 r

1

2

(b) Wendland kernel 2

Quintic First derivative Seconnd derivative

100

Improved Gaussian First derivative Seconnd derivative

1.5 1

50

W/

W/

0.5 0

0 -0.5

-50

-1 -1.5

-100 -3

-2

-1

0 r

1

(c) Quintic spline kernel

2

3

-2 -3

-2

-1

0 r

1

2

(d) Improved Gaussian kernel

Fig. 1. The kernel functions and their derivatives.

3

244

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

0.080 Cubic Wendland Quintic Improved Gaussian

0.060 0.040 0.020 0.000 -0.020 0

0.1

0.2 0.3 d/dx

0.4

0.5 Fig. 4. The initial pressure ﬁeld.

Fig. 2. Curves of the function approximation errors.

and their derivatives of the ﬁrst two orders are smooth enough, but the second order derivatives of the cubic spline kernel and the Wendland kernel have inﬂection points. The effect of these four kernel functions on accuracy and computational results will be further investigated later. Taking the two-dimensional function for example, the function and its derivatives are f (x, y) = x2 + y2

(4)

fx (x, y) = 2x, fy (x, y) = 2y

The two-dimensional problem domain is assumed to be x ∈ [0, 2] and y ∈ [0, 2], the uniform particles spacing dx = 0.1 and dy = 0.1, and the smoothing length h = 1.23dx. The distance between the probe and the center of the domain is noted as d. Comparisons of the particle approximation errors of the above four kernel functions are shown in Figs. 2 and 3, where =<

f (x, y)

>

vx0 = −A0 x

(6)

vy0 = A0 y and the initial pressure ﬁeld

(7)

(5)

−f (x, y)

Because the derivatives with respect to x and y have the same results, only the approximation errors of the function’s derivative with respect to x are presented. As seen from the ﬁgures, the approximation errors of the functions and their derivatives are almost the same for the improved Gaussian kernel and the quintic spline kernel. The accuracy of the cubic spline kernel derivative attenuates more seriously. A large error occurs in the approximation by using Wendland kernel, which accords with the conclusion of [50]. It pointed out in [50] that the Wendland kernel

1.000 Cubic Wendland Quintic Improved Gaussian

0.500

0.000

-0.500

-1.000 0

P0 (x, y) = 0.5A0 [R2 − (x2 + y2 )]

ε =< f (x, y) > −f (x, y) ε

function may cause large errors with less neighbor particles, and about 21 particles are used in this simulation. When solving the hydrodynamics problems, partial differential equations about density, velocity and other variables need to be solved; therefore the derivative of a kernel function has great signiﬁcance. The approximation errors of the function’s derivative are less than 0.8% for the improved Gaussian kernel and the quintic spline kernel with uniform particle distribution. In the following contents, a ﬂow example is taken to investigate the computational stability of the four commonly used kernel functions. There is a circle basin of radius R, with the initial velocity ﬁeld

0.1

0.2 0.3 d/dx

0.4

0.5

Fig. 3. Curves of the function’s derivative approximation errors.

The initial pressure distribution is shown in Fig. 4. The particles are distributed uniformly at initial time. The motion is in two-dimensional plane without any external force. Because of the symmetry, only the upper part of the pressure ﬁeld is presented. The pressure distributions for different kernel functions are shown in Fig. 5. As seen from Fig. 5, the result of cubic kernel function is unstable and the instability occurs. The result of WC2 shows big errors, similar to the previous case. As for the initial velocity amplitude A0 and 2A0 , more stable results are obtained from the improved Gaussian kernel function and the quintic spline kernel function. The results of these two kernel functions are almost the same when the circle basin reaches the same displacement. They also show that the Gaussian kernel function and the quintic spline kernel function are not sensitive to the disorder distribution of particles. Comparing the above results, the improved Gaussian kernel and the quintic spline kernel have better accuracy and more stable results. The approximation errors of the derivatives are only 0.8% for the ﬁxed particle case, and the two kernel functions are not sensitive to the disorder particle distribution. The compact support domain of the improved Gaussian kernel and the quintic spline kernel is also the same. Besides, it is convenient and efﬁcient to obtain the derivative of Gaussian kernel function as ∂W ∂xid

= −2

xid − xjd h2

Wij

(8)

Therefore the improved Gaussian kernel is selected in this paper.

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

245

Fig. 5. Pressure ﬁeld of different velocities with different kernel functions.

2.3. Boundary treatment When dealing with hydrodynamics problems, several types of boundaries always encountered, such as the solid wall boundary and the free surface. Because of the existence of boundaries, the support domain of particles near boundaries is truncated, which will introduce calculation errors and affect the accuracy of numerical simulations. As for the free surface, most of SPH literatures state that there is no need for a special treatment. At ﬁrst, the SPH method is a Lagrangian method using Lagrangian particle dispersion, so the kinematic conditions of free surface can be considered as automatically satisﬁed. As for the dynamic boundary conditions of free surface, although the pressure at the free surface is not strictly zero in the actual calculations, it is guaranteed that the governing equation converges to order h [49] when the smoothing length approaches zero. How to deal with the solid wall boundary in SPH method is a difﬁcult and focused issue. Many methods have been put forward. In summary, there are two main sorts of methods. One lays multilayer particles outside the boundary to ensure that the support domain of ﬂuid particles near the boundary is not attenuated by the wall particles; the other lays a single layer of particles at the boundary, so that repulsive force is generated to avoid particles penetration when the ﬂuid particles approach the boundary. There are also many combined models for boundary particles and ﬂuid particles utilizing the above two ideas [31,52–57]. But so far, there is no better way to unify the boundary treatment. In the repulsive boundary model, the non-physical disturbance always exists in the ﬂow ﬁeld. In this paper, the new formulation for boundary condition by Adami et al. [54] is selected to calculate the pressure of particles on the solid wall boundary. At the initial time, the boundary particles are well arranged once and for all, and the relative position is not changed during the calculation. It could be drawn from the results by Adami et al. [54] that, this treatment is not only capable of simulating complex boundary, but also easy to be extended to three-dimensional problems. Some relatively stable results are obtained through this boundary method. The particle distribution of the solid wall boundary is shown in Fig. 6. Four layers of particles are laid to model the solid boundary to avoid numerical errors caused from boundary truncation. The physical quantities of ﬂuid particles, such as density and velocity, are obtained by solving Eq. (1). The support domain of ﬂuid particles includes the interacting ﬂuid particles and the boundary particles. The pressure of boundary particles is obtained by interpolation from the adjacent ﬂuid particles but not by solving the continuity equation and the equation of state. The density of boundary particles can be obtained by solving the equation of state (2). For the test cases in this paper, the liquid tanks are under forced motion. The velocity of dummy boundary particles is directly

Fig. 6. Boundary particle distribution.

interpolated from the given tank velocity. According to the force balance between the ﬂuid particle and boundary particle, one can get: af · r bf = ab · r bf

(9)

where, a and r denote the acceleration and position; besides, rbf = rb − rf , the subscripts b and f represent the solid wall boundary particles and ﬂuid particles respectively. According to the momentum equation, Eq. (9) can be transformed to

∇ P · r bf = f (g − ab ) · r bf

(10)

For a pair of interacting particles, ∇ P · r bf is equal to the pressure increment along the centerline of boundary particle and ﬂuid particle. Thus, the pressure of boundary particle can be drawn from Pb − Pf = f (g − ab ) · r bf

(11)

by the SPH approximation of Eq. (11), the pressure of boundary particles is expressed as: Pb =

(Pf + f (g − ab ) · r bf )Wbf

(12)

f

Since the boundary particle information can only be obtained through the ﬂuid particles, the normalization of Eq. (12) is always carried out to overcome the defects of boundary, and the ﬁnal expression of the boundary particle pressure is

Pb =

f

(Pf + f (g − ab ) · r bf )Wbf

f

Wbf

(13)

After obtaining the pressure of boundary particles through Eq. (13), the particle density can be got by solving the Eq. (2): b = 0

P

b

B

1/

+1

(14)

The pressure formula of boundary particles considers not only the inﬂuence of the interacting ﬂuid particles, but also the acceleration of the boundary particles themselves. Also, it can be applied to the moving solid boundary.

246

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

Fig. 7. Pressure distribution of a dam break.

Fig. 7 shows the pressure distribution of a benchmark of dam break with the boundary used in this paper. As plotted in the ﬁgure, the pressure of boundary particles obtained from the interacting ﬂuid particles shows a good continuity and no disturbance. In this paper, the predictor-corrector method is chosen as the time integration algorithm [58], and the program is paralleled [60]. In order to improve the continuity, the continuity correction of density ﬁeld is carried out using the moving least squares (MLS) method [17] every 20 steps. The MLS correction satisﬁes the zeroth order and ﬁrst order continuity conditions throughout the discrete domain, and has ﬁrst-order completeness [59]. Although the MLS correction makes the calculation slightly expensive, it can improve the accuracy and stability of the calculation. The ﬂowchart of the SPH program is shown in Fig. 8, where tmax is the maximum time of numerical simulation.

3. Results and discussions The researches of liquid sloshing in tanks are always focused on the impact pressure on the tank, the inﬂuences of internal structures, the coupling problems of multi-degrees of freedom and so on. Therefore, the following three models are chosen to validate the numerical methods in this paper. The distribution of the particles in the computational models is uniform. In the two-dimensional models, the particle spacing is dx = 0.005 m and the smoothing length is h = 1.23dx. A constant time step dt = 4 × 10−5 s is employed. In the three-dimensional models, the particle spacing is dx = 0.01 m and the smoothing length is h = dx. A constant time step dt = 8 × 10−5 s is adopted. The tanks are partially ﬁlled with water and the effects of air and structural deformation are not considered in the simulation. For a given rectangular tank partially ﬁlled with water, its natural frequency is ωn2 = g(n/L) tanh(n(D/L), where, Land D represents the width of the tank and the depth of the liquid; n is the modulus. ω1 (n = 1) is the ﬁrst order natural frequency, namely the lowest natural frequency. The liquid sloshing violently at the lowest natural frequency and the maximum impact pressure also occurs around this frequency. 3.1. Rolling motion of a rectangular tank In this case, the rolling motion of a rectangular tank is studied, as shown in Fig. 9. The tank model is the same as Akyildiz’s [6] experimental model: the rectangular tank has a length of L = 0.92 m, a height of H = 0.62 m, a liquid depth of D = 0.25 H, and thus the lowest natural frequency is ω1 = 4.025 rad/s. The external excitation is = 0 sin(ωr t), where 0 is the angular displacement and ωr is the circular frequency of the rolling motion. The rolling axis is located

Fig. 8. Flowchart of the SPH program.

Fig. 9. Model of a rolling tank.

Table 2 Calculation cases.

Case 1 Case 2 Case 3

0 (◦ )

ωr (rad/s)

4 4 4

2 4 6

Case 4 Case 5 Case 6

0 (◦ )

ωr (rad/s)

8 8 8

2 4 6

Case 7 Case 8 Case 9

0 (◦ )

ωr (rad/s)

12 12 12

2 4 6

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

247

Fig. 10. Pressure distribution with 0 = 8◦ at ωr = 2.0 rad/s.

Fig. 11. Pressure distribution with 0 = 8◦ at ωr = 2.0 rad/s, ωr = 4.0 rad/s, ωr = 6.0 rad/s.

2

2.5

Experiment SPH D/dx=15 D/dx=40 D/dx=60

2

1.5

P/kPa

P/kPa

1.5

Experiment SPH D/dx=15 D/dx=40 D/dx=60

1

1

0.5 0.5

0 0

5

10

15

20

0 0

5

10

(a)

0

4

15

t/s

t/s

(b) Fig. 12. Comparison of the probe pressure.

0

8

20

248

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

6 =4 =8 =12

5

P/kPa

4 3 2 1 0

2

3

4

5

6

/ rad s-1 r Fig. 13. Comparison of the maximum pressure.

Fig. 14. Model of a tank with a bafﬂe.

Fig. 15. Typical velocity for model (3).

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

249

Fig. 16. Typical velocity for models (1), (2) and (4).

at the central axis of the model; the measuring probe is located at the wall, 0.06 m away from the bottom. The calculation cases are shown in Table 2. The case where 0 = 8◦ is taken as an example to analyze the sloshing wave motion in the tank. The law of sloshing motion in the ﬁrst period is shown in Fig. 10, with the rolling amplitude being 0 = 8◦ . With the rightward motion of the tank from the middle position, the liquid impacts on the tank and the wall pressure increases; when the angle increases to the maximum, the velocity of the tank decreases to zero and the tank starts the reverse movement, while the liquid is still in the rightward motion, which makes the sloshing wave impact on the leftward moving wall and results in the leftward movement of the wave. The wave ﬁrstly arrives at the left wall and then starts rightward movement after the collision with the left wall, while the tank continues its leftward rolling until the maximum amplitude and starts the reverse movement. In the initial period of the rolling tank, there is a certain phase differences

between the motion of tank and the sloshing wave; however, with the forced rolling motion of the tank, the motion of sloshing wave catches up with the motion of the tank gradually. Fig. 11 shows the pressure distribution at different excitation frequencies. When the amplitude is 0 = 8◦ and ωr =4 rad/s, the liquid sloshing is the most violent. Fig. 12 shows the comparison of the pressure on measuring probe from the experimental results and the numerical results at different excitation frequencies ( 0 = 4◦ and 0 = 8◦ ) within 20 s. As seen from the ﬁgure, the numerical results of the two models are in good agreements with the experiment results. Two additional SPH results are also shown in the ﬁgure with D/dx being 15 and 60 while h/dx kept the same. The convergence of the SPH results is shown in the ﬁgure. The numerical method in this paper shows stable results. With the increase of the sloshing excitation amplitude, the impact pressure on the measuring probe increases, and the pulse time of impact pressure is reduced. Fig. 13 shows the comparison between

250

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

0.03 Experiment SPH result

0.02

h/m

0.01 0 -0.01 -0.02 0

5

10

15

20

25

t/s Fig. 17. Comparison of the wave height for model (3).

Fig. 18. Comparison of wave height for P2 .

the average peak pressures on the measuring probe at nine different calculation cases from10 s to 20 s. As seen from the ﬁgure, as the rolling angle increases, the peak of the impact pressure increases signiﬁcantly; the impact pressure does not increase with the sloshing frequency, but increases signiﬁcantly around the lowest natural frequency, in which case the pressure on the measuring probe can be as much as twice of the pressure values at other frequencies. 3.2. Surging motion of a rectangular tank with a bafﬂe The model of the sloshing tank is shown in Fig. 14. The tank has a length of L = 0.60 m and a height of H = 0.60 m. The depth of the

Fig. 19. Comparison of pressure for P1 .

ﬁlled ﬂuid is D = 0.30 m. A vertical bafﬂe is mounted at the middle of the tank bottom, whose height is d. A measuring probe P1 is located on the wall in the still free surface, and another measuring probe P2 is located at 0.04 m away from the wall. The surging motion of a rectangular tank is studied in this case with an excitation of y = Asin(ωy t), an amplitude of A = 0.0024 m, and a frequency of ωy = ω1 . The height of the bafﬂe is d/D = 0 for model (1); d/D = 0.30 for model (2), d/D = 0.50 for model (3), and d/D = 0.70 for model (4). The typical velocity ﬁled for model (3) d/D = 0.50 and the velocity vector near the bafﬂe are shown in Fig. 15. The vortices evolution and vortex shedding are also plotted. Because the bafﬂe width is very small compared with the tank width, it is difﬁcult to capture the detailed ﬂow phenomenon near the bafﬂe. As seen from the ﬁgure, the SPH method has successfully captured the vortices evolution and vortex shedding near the bafﬂe. When the tank moves from right to left, it drives the sloshing waves to move from right to left, which will generate counterclockwise vortex on the left of the bafﬂe tip, as shown in Fig. 15(a). On the contrary, when the tank moves from left to right, it will generate clockwise vortex on right of the bafﬂe tip, as presented in Fig. 15(d). With the motion of sloshing waves, the vortex will move upward, which can be seen in Fig. 15(b) and (e), and gradually disappear, as plotted in Fig. 15(c) and (f). Fig. 16 shows the typical velocity for model (1), model (2) and model (3). Comparing the four numerical models, the maximum velocity of model (1) without bafﬂe is on the free surface, but for the models with a bafﬂe, the maximum velocity is mainly at the bafﬂe tip. Because of the bafﬂe, vortex near the bafﬂe tip is generated, which will cause energy dissipative, and thereby reduce the

Fig. 20. Pressure distribution for model (3).

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

251

Fig. 21. Model of three-dimensional tanks.

swaying excitation is x = Asin(ωx t)sin45◦ and the surging excitation is y = Asin(ωx t)cos45◦ . Fig. 22 shows the comparison between wave height for model (1) and experimental results [1] in 18 s. The numerical result is in good agreements with the experiment trends. In the simulation the external excitation is inserted at the starting time t = 0 s. There are some differences in the amplitude, which is reasonable since the liquid sloshing in tank is a complex ﬂuid motion characterized by nonlinearity and strong randomness, and the tank motion is coupled with swaying and surging. In Fig. 23, the comparison of the wave height on the probe for model (1) and model

Experiment SPH result

0.04

0.02

h/m

kinetic energy of the free surface and effectively reduce the height of sloshing wave. Comparisons of the wave height for numerical and experimental results [9] at 25 s are shown in Fig. 17. The SPH numerical results are in good agreements with the experimental results. The excitation frequency is the same as the lowest natural frequency of the tank and thus the liquid sloshing is violent, but quickly attenuated to a stable value. The comparison of wave height on P2 and impact pressure value on P1 between 10 s and 20 s for the four models are presented in Figs. 18 and 19. As shown in Fig. 18, the bafﬂe effectively reduces the height of the sloshing wave; the maximum wave height for model (4) is only 6.5% of that for model (1). With the increasing of bafﬂe height, the wave height on the probe decreases. It is the same for the pressure on P1 in Fig. 19, and the maximum wave height for model (4) is about 15% of that for model (1). From the two ﬁgures, it can be found when d/D ≤ 0.50, the values on the probes are similar, but when d/D > 0.50, the pressure on P1 and the wave height on the probe P2 are reduced by more than a half. The typical pressure distributions are given in Fig. 20, because of the existence of bafﬂe, the pressure ﬁeld of ﬂuid is discontinuous near the bafﬂe, and the change of pressure ﬁeld is not signiﬁcant. From the above results it can see that increasing the bafﬂe height to over a certain percentages of the liquid height is equivalent to reducing the tank width, and the lowest natural frequency of the tank will be affected. From Section 3.1, one knows that when the sloshing frequency is far away from the lowest natural frequency of the tank, the impact pressure is signiﬁcantly reduced. Therefore, the bafﬂe effectively reduces the impact of the liquid on the tank.

0

-0.02

-0.04 0

5

10

15

t/s Fig. 22. Comparison of wave height for model (1).

3.3. Coupling motion of a three-dimensional sloshing tank

Without baffle With baffle

0.04

0.02

h/m

A model can be simpliﬁed into a two-dimensional model for sloshing with a single degree of freedom, but for the coupling motion of multi-degrees of freedom, it is essential to establish three-dimensional models. In this case, the coupling motion of swaying and surging of a three-dimensional sloshing tank is studied. The tank has a size of L × B × H = 0.59 m × 0.59 m × 0.59 m, the depth of ﬂuid is D = 0.508H. Model (1) is without the bafﬂe but model (2) has a vertical bafﬂe mounted on the middle of the tank bottom, whose height is d = 0.5D. There are three measuring probes in the numerical models, as plotted in Fig. 21: the probe P1 is located at the middle of the tank bottom and 0.04 m away from the wall; P2 is located on the wall in the still free surface along the excitation direction; P3 is at the center of the wall which is parallel to the bafﬂe and 0.15 m away from the bottom. The model of the threedimensional tanks is shown in Fig. 21. The external excitation has an amplitude of A = 0.0078L, and a frequency of ωx = ωy = 1.1ω1 ; the

0

-0.02

-0.04 0

5

10

15

t/s Fig. 23. Comparison of wave height for model (1) and model (2).

252

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

Fig. 24. Height of the free surface of model (1).

Fig. 25. Height of the free surface of model (2).

(2) at 18 s is shown. The vertical bafﬂe in the middle of the tank bottom is effective to reduce the height of the sloshing wave. The free surface deformation for the above two models are shown in Figs. 24 and 25 respectively, where the thicker black line represents the outline of the liquid tank and the thinner one is the

initial height of the still water. For the model without bafﬂe, due to the diagonal excitation, the sloshing wave is also along the diagonal direction. In Fig. 25, the bafﬂe affects not only the height but also the direction of the sloshing wave. The difference in wave height for the two models is about 0.03 m, and the magnitude of the liquid

2

0.8 Without baffle With baffle

1.8

P / kPa

P / kPa

0.6

0.4

Without baffle With baffle

1.6 1.4

0.2 1.2 0 0

5

10

t/s

(a) Pressure on P2

15

1 0

5

10

15

t/s

(b) Pressure on P3

Fig. 26. Comparisons of the pressure for model (1) and model (2).

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

sloshing reduces about 50% for d/D = 0.5. The direction of the sloshing wave is also changed; the wave has a tendency to transform from diagonal to planar, which moves from one side of the tank wall to the opposite side. The three-dimensional sloshing problems also show good results. The pressure curves at probes 2 and 3 for the above models are shown in Fig. 26. It is obvious that the bafﬂe has some attenuation effects on the impact of inner liquid. However, compared with that in the motions of a single degree of freedom, the effects of the bafﬂe on the impact pressure reduce greatly in the coupling motions. Moreover, the pressure at probe 3 becomes more complicated in the coupling motion. The bafﬂe is more effective for reducing the sloshing wave height than the impact pressure. 4. Conclusions Firstly, several core issues of SPH method, the kernel function and the boundary treatment, are investigated: in term of the accuracy and stability of the kernel function, the Gaussian kernel function is proved to be suitable for the present hydrodynamics problems; in order to solve the truncation of the support domain of the kernel function near boundary, the dummy particles boundary is introduced and meanwhile the boundary motion is considered. Afterwards, the paralleled SPH codes is programed to study the liquid sloshing in both two-dimensional and three-dimensional tanks of single and multi-degrees of freedom. For the rolling motion of the tank, with the increase of the excitation amplitude, the impact pressure on the tank is increased, and the peak of impact pressure doubles when excitation frequency is close to the lowest natural frequency of the tank. For the tank with a bafﬂe in surging motion, the bafﬂe can effectively reduce the height of sloshing wave by more than 90% in two-dimensional tank under a single excitation; when the ratio of the bafﬂe height to the still liquid height is less than 0.5, the variation of the lowest natural frequency is not obvious; however, when the ratio is more than 0.5, the liquid sloshing in tank reduces signiﬁcantly, this is equivalent to reducing the tank width and will increase the lowest natural frequency of the tank. For the three-dimensional sloshing tank under a coupled excitation, the bafﬂe can change the direction of sloshing wave and at the same time the magnitude of the liquid sloshing is reduced by about 50%. But the reduction of the sloshing impact pressure is not obvious. All the numerical results in this paper are in good agreements with the existing experimental results. To sum up, the reliability and stability of the parallel program are veriﬁed which lays a foundation for further studies for the sloshing coupled with ship motions. Acknowledgements This work is supported by the National Program for Special Support of Eminent Professionals, the Excellent Young Scientists Fund (51222904), the Excellent Youth Foundation of Heilongjiang Province (JC201207), and Lloyd’s Register Foundation, which helps to protect life and property by supporting engineering-related education, public engagement and the application of research. References [1] Faltinsen OM, Rognebakke OF, Timokha AN. Resonance three-dimensional nonlinear sloshing in a square-base basin. J Fluid Mech 2003;487:1–42. [2] Faltinsen OM, Rognebakke OF, Timokha AN. Resonant three-dimensional nonlinear sloshing in a square-base basin. Part 2. Effect of higher modes. J Fluid Mech 2005;523:199–218. [3] Faltinsen OM, Rognebakke OF, Timokha AN. Resonant three-dimensional nonlinear sloshing in a square-base basin. Part 3. Base ratio perturbations. J Fluid Mech 2006;551:93–116. [4] Huang S, Duan WY, Ma QW. An approximation to energy dissipation in time domain simulation of sloshing waves based on linear potential theory. China Ocean Eng 2011;2:189–200.

253

[5] Ikeda T, Nakagawa N. Non-linear vibrations of a structure caused by water sloshing in a rectangular tank. J Sound Vib 1997;201(1):23–41. [6] Akyildiz H, Ünal E. Experimental investigation of pressure distribution on a rectangular tank due to the liquid sloshing. Ocean Eng 2005;32(11):1503–16. [7] Pistani F, Thiagarajan K. Experimental measurements and data analysis of the impact pressures in a sloshing experiment. Ocean Eng 2012;52:60–74. [8] Molin B, Remy F. Experimental and numerical study of the sloshing motion in a rectangular tank with a perforated screen. J Fluid Struct 2013;43:463–80. [9] Wu CH, Chen BF, Hung TK. Hydrodynamic forces induced by transient sloshing in a 3D rectangular tank due to oblique horizontal excitation. Comput Math Appl 2013;65(8):1163–86. [10] Nasar T, Sannasiraj SA, Sundar V. Experimental study of liquid sloshing dynamics in a barge carrying tank. Fluid Dyn Res 2008;40(6): 427–58. [11] Liu DM, Lin PZ. A numerical study of three-dimensional liquid sloshing in tanks. J Comput Phys 2008;227:3921–39. [12] Nagashima T. Sloshing analysis of a liquid storage container using level set X-FEM. Commun Numer Methods Eng 2009;25(4):357–79. [13] Lee D, Kim M, Kwon S, Lee Y. A parametric sensitivity study on LNG tank sloshing loads by numerical simulations. Ocean Eng 2007;34:3–9. [14] Wu CH, Faltinsen OM, Chen BF. Numerical study of sloshing liquid in tanks with bafﬂes by time-independent ﬁnite difference and ﬁctitious cell method. Comput Fluids 2012;63:9–26. [15] Colagrossi A, Antuono M, Le Touzé D. Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model. Phys Rev E 2009;79(5):056701. [16] Monaghan JJ. Simulating free surface ﬂows with SPH. J Comput Phys 1994;110(2):399–406. [17] Colagrossi A, Landrini M. Numerical simulation of interfacial ﬂows by smoothed particle hydrodynamics. J Comput Phys 2003;191(2):448–75. [18] Molteni D, Colagrossi A. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Comput Phys Commun 2009;180(6):861–72. [19] Antuono M, Colagrossi A, Marrone S, Molteni D. Free-surface ﬂows solved by means of SPH schemes with numerical diffusive terms. Comput Phys Commun 2010;181(3):532–49. [20] Antuono M, Colagrossi A, Marrone S. Numerical diffusive terms in weaklycompressible SPH schemes. Comput Phys Commun 2012;183(12):2570–80. [21] Raﬁee A, Cummins S, Rudman M, Thiagarajan K. Comparative study on the accuracy and stability of SPH schemes in simulating energetic free-surface ﬂows. Eur J Mech B Fluids 2012;36:1–16. [22] Sirotkin FV, Yoh JJ. A smoothed particle hydrodynamics method with approximate Riemann solvers for simulation of strong explosions. Comput Fluids 2013;88:418–29. [23] Ferrari A, Dumbser M, Toro EF, Armanini A. A new 3D parallel SPH scheme for free surface ﬂows. Comput Fluids 2009;38(6):1203–17. [24] Roubtsova V, Kahawita R. The SPH technique applied to free surface ﬂows. Comput Fluids 2006;35(10):1359–71. [25] Hu XY, Adams NA. An incompressible multi-phase SPH method. J Comput Phys 2007;227(1):264–78. [26] Ming FR, Zhang AM, Cao XY. A robust shell element in meshfree SPH method. Acta Mech Sin 2013;29(2):241–55. [27] Zhang AM, Ming FR, Wang SP. Coupled SPHS–BEM method for transient ﬂuid–structure interaction and applications in underwater impacts. Appl Ocean Res 2013;43:223–33. [28] Zheng X, Duan WY, Ma QW. A new scheme for identifying free surface particles in improved SPH. Sci China Phys Mech Astron 2012;55(8):1454–63. [29] Omidvar P, Stansby PK, Rogers BD. SPH for 3D ﬂoating bodies using variable mass particle distribution. Int J Numer Methods Fluids 2013;72(4):427–52. [30] Delorme L, Colagrossi A, Souto-Iglesias A, Botia-Vera E. A set of canonical problems in sloshing, Part I: Pressure ﬁeld in forced roll – comparison between experimental results and SPH. Ocean Eng 2009;36(2):168–78. [31] Shao JR, Li HQ, Liu GR, Liu MB. An improved SPH method for modeling liquid sloshing dynamics. Comput Struct 2012;100:18–26. [32] Marsh A, Prakash M, Semercigil E, Turan ÖF. A study of sloshing absorber geometry for structural control with SPH. J Fluid Struct 2011;27(8):1165–81. [33] Colagrossi A, Lugni C, Brocchini M. A study of violent sloshing wave impacts using an improved SPH method. J Hydraul Res 2010;48(S1):94–104. [34] Khayyer A, Gotoh H. Wave impact pressure calculations by improved SPH methods. Int J Offshore Polar Eng 2009;19(4):300–7. [35] Gotoh H, Khayyer A, Ikari H, Arikawa T, Shimosako K. On enhancement of Incompressible SPH method for simulation of violent sloshing ﬂows. Appl Ocean Res 2014;46:104–15. [36] Vorobyev A, Kriventsev V, Maschek W. Simulation of central sloshing experiments with smoothed particle hydrodynamics (SPH) method. Nucl Eng Des 2011;241(8):3086–96. [37] Raﬁee A, Pistani F, Thiagarajan K. Study of liquid sloshing: numerical and experimental approach. Comput Mech 2011;47(1):65–75. [38] Shen L, Vassalos D. Applications of 3D parallel SPH for sloshing and ﬂooding. Contemporary Ideas on Ship Stability and Capsizing in Waves. Netherlands: Springer; 2011. p. 709–21. [39] Bai ZG, Zhao J, Zhang W, Wang WL. 3D SPH numerical investigation for the sloshing impact in LNG tank. In: 32nd International Conference on Ocean, Offshore and Arctic Engineering OMAE. 2013. [40] Fey M, Jeltsch R. Hyperbolic problems: theory, numerics, applications. Int Ser Numer Math 1999;130:897–8.

254

X.Y. Cao et al. / Applied Ocean Research 47 (2014) 241–254

[41] Liu GR, Liu MB. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientiﬁc; 2003. [42] Allahdadi FA, Carney TC, Hipp JR, Petschek AG. High strain Lagrangian hydrodynamics: a three dimensional SPH code for dynamic material response. Phillips Lab Kirtland AFB NM; 1993. [43] Monaghan JJ. Smoothed particle hydrodynamics. Rep Prog Phys 2005;68(8):1703–59. [44] Cole RH, Weller R. Underwater explosions. New Jersey: Princeton University Press; 1948. [45] Schonberg IJ. Contributions to the problem of approximation of equidistant data by analytic functions. Quart Appl Math 1946;4(2):45–99. [46] Colagrossi A, Souto-Iglesias A, Antuono M, Marrone S. Smoothed-particlehydrodynamics modeling of dissipation mechanisms in gravity waves. Phys Rev E 2013;87(2):023302. [47] Liu X, Xu H, Shao S, Lin P. An improved incompressible SPH model for simulation of wave–structure interaction. Comput Fluids 2013;71:113–23. [48] Asai M, Aly AM, Sonoda Y, Sakai Y. A stabilized incompressible SPH method by relaxing the density invariance condition. J Appl Math 2012;2012:139583. [49] Marrone S. Enhanced SPH modeling of free-surface ﬂows with large deformations. 2012. [50] Dehnen W, Aly H. Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Mon Not R Astron Soc 2012;425(2):1068–82.

[51] Yang XF, Liu MB. An improvement for stress instability in smoothed particle hydrodynamics. Acta Phys Sin 2012;61(22):224701. [52] Marrone S, Antuono M, Colagrossi A, Colicchio G, Le Touzé D, Graziani G. ␦-SPH model for simulating violent impact ﬂows. Comput Methods Appl Mech Eng 2011;200(13):1526–42. [53] Oger G, Doring M, Alessandrini B, Ferrant P. Two-dimensional SPH simulations of wedge water entries. J Comput Phys 2006;213(2):803–22. [54] Adami S, Hu XY, Adams NA. A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys 2012;705(231): 7–7075. [55] De Leffe M, Le Touzé D, Alessandrini B. Normal ﬂux method at the boundary for SPH. In: Proceedings of the 4th SPHERIC Workshop. 2009. p. 150–7. [56] Rogers BD, Dalrymple RA. SPH modeling of tsunami waves. Coastal Ocean Eng 2008;10:75–100. [57] Monaghan JJ, Kajtar JB. SPH particle boundary forces for arbitrary boundaries. Comput Phys Commun 2009;180(10):1811–20. [58] Monaghan JJ. On the problem of penetration in particle methods. J Comput Phys 1989;82:1–15. [59] Belytschko T, Krongauz Y, Dolbow J, Gerlach C. On the completeness of meshfree particle methods. Int J Numer Methods Eng 1998;43:785–819. [60] Zhang AM, Cao XY, Ming FR, Zhang ZF. Investigation on a damaged ship model sinking into water based on three-dimensional SPH method. Appl Ocean Res 2013;42:24–31.