MING Ping-jian

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China, E-mail: pingjianming@hrbeu.edu.cn

DUAN Wen-yang

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China

NUMERICAL SIMULATION OF SLOSHING IN RECTANGULAR TANK WITH VOF BASED ON UNSTRUCTURED GRIDS*

MING Ping-jian

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China, E-mail: pingjianming@hrbeu.edu.cn

DUAN Wen-yang

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China

(Received June 15, 2010, Revised August 30, 2010)

A new method for sloshing simulation in a sway tank is present, in which the two phase interface is treated as a physical discontinuity, which can be captured by a well-designed high order scheme. Based on Normalized Variable Diagram (NVD), a high order discretization scheme with unstructured grids is realized, together with a numerical method for free surface flow with a fixed grid. This method is implemented in an in-house code General Transport Equation Analyzer ( GTEA ) which is an unstructured grids finite volume solver. The present method is first validated by available analytical solutions. A simulation for a 2-D rectangular tank at different excitation frequencies of the sway is carried out. A comparison with experimental data in literature and results obtained by commercial software CFX shows that the sloshing load on the monitor points agrees well with the experimental data, with the same grids, and the present method gives better results on the secondary peak. It is shown that the present method can simulate the free surface overturning and breakup phenomena.

sloshing, numerical simulation, unstructured grid, Volume Of Fluid (VOF)

1. Introduction

Sloshing occurs in partially filled tank subjected to external excitations or with fluid of changing states. The larger the volume of a ship, the more possible its exposure to excessive sloshing loads during the ship operation life. In recent years, the demand for Liquefied Natural Gas (LNG) is increasing for the environmental and other reasons. It is an urgent task to develop new generation LNG ships with super-large capacity, from the current 140 000 m3or so to 250 000 m3. For such a large LNG ship, an important prerequisite for its safe operation is an accurate prediction of the sloshing loads, therefore the slamming load calculation for tank liquid sloshing attracts a worldwide attention. As the fluid sloshing is in essence a non-linear movement, there is no theoretical solution so far. The experimental and numerical simulations are two main methods at present. Sloshing problems have received extensive attentions over the past few years, and a lot of research methods have been proposed. Experimental studies provide verifications for the numerical and theoretical methods[1-3]. Sloshing problems are solved based on potential theory, analytical method[4,5], the boundary element method[6]and the finite element method[7,8]. In many situations, the fluid viscosity effect cannot be neglected and the Navier-Stokes equation has to be used to describe sloshing. Hirt and Nicols firstly proposed the Volume Of Fluid (VOF) method, in which the free surface is constructed by segments parallel to the coordinates, which was applied widely subsequently. It is called Simple Linear Interface Construction (SLIC). Youngs improved the SLIC method and put forward thePiecewise Linear Interface Construction (PLIC), with good results[9]. Yang simulated the interaction between extreme wave and freely-floating structure with VOF method based on unstructured grid[10]. Ai proposed an unstructured grid method to simulate non-hydrostatic free surface flow[11]. But neither method can guarantee that the free surface on either side of the cell face calculated would completely overlap. So the calculation result may be difficult to converge, or even to diverge. A number of improved methods were proposed, with good results[12]. A new SPH method was adapted to simulate the water wave[13]. The advantages and disadvantages of various numerical methods are reviewed[14,15].

With the VOF mentioned above, the free surface is constructed by geometric method. However, the reconstruction of the surface takes up a large amount of computation resources, mainly based on structured grids. This article proposes a new concept that the free surface is treated as a physical discontinuity, which can be captured by a well designed high-order scheme the same as the shock captured scheme on unstructured grids. Therefore the sloshing process can be simulated on fixed unstructured grids.

This article is organized as follows. In Section 2, the mathematical models and numerical methods are presented. In Section 3, computation models and the analytical solution are given. Section 4 shows the results and contains a discussion. Conclusions are in the final Section 5.

2. Mathematical models and numerical methods

2.1 Mathematical models

The Navier-Stockes equations are used in this article as the governing equations for sloshing. A homogeneous model is adopted in which the two phases share the same pressure and velocity. The gas-liquid interface is looked as a physical discontinuity where the properties change sharply. Hence, the control equation is the same as a single-phase flow.

Physical parameters in the two-phase flow can be calculated as follows

where α is defined as

The equation of state is needed to determine the relationship between pressure and density and to close the control equations for compressible flow, such as the ideal gas EOS as

For incompressible flow, ρ is constant, Eq.(1) can be written in the following form

In all calculations without special statement in this article, it is assumed that both phases are incompressible.

For the tank sway process, the non-inertial coordinates (XOY) show in Fig.1 are used, with the origin at tank center. The movement of the tanks in the XOY coordinate system is described as

The acceleration is

The movement of the coordinate system XOY at O point is

The acceleration

Therefore, the volume force is

2.2 Numerical method

A second order Crank-Nicolson scheme is used for time integral discretization and the unstructured finite volume method for spatial terms. A Semi-Implicit Method for Pressure Link Equations ( SIMPLEs ) algorithm is adopted for the pressure velocity coupling. It is a prediction-correction two stage method.

(1) Momentum prediction

Assuming or knowing the pressure distribution, we solve Eq.(2) to obtain the velocity distribution. At this time, the velocity satisfies the momentum conservation equation, but may not satisfy the mass conservation equation. The velocity field needs to be corrected according to the mass conservation equation.

(2) Pressure-correction

Unfortunately, there is no control equation for pressure. The pressure velocity relation should be derived from the moment equation and then a control equation for pressure correction can be obtained based on the mass conservation. We can correct the velocity and update the pressure field at the same time.

The above procedure is used in many studies and in almost all popular CFD commercial softwares. The special step for sloshing is the solution of the volume fraction equation, which will be discussed below in details.

2.3 NVD and high-order scheme

For a collocated grid finite volume method, all variables are referred to the cell center. For the convection term, variable values on the cell face are required and there are many different interpolation method based discretization schemes.

For clearness, 1-D case is used to explain the Normalized Variable Diagram ( NVD) procedure[16]. Based on the NVD, several TVD schemes can be designed[17]. As shown in Fig.2, the cell face value

fα

is under consideration. According to the flow direction of the cell face, we can define its upstream and downstream cell valuesDαandAα, and the next upstream valueUα, and then dimensionless variables can be defined accordingly,

Therefore, the dimensionless variables of the upstream cell take the form

And the cell interface dimensionless variables take the form

Different discretization forms can be expressed by NVD as shown in Fig.3[18].

Different discretization forms can be expressed by the following relationship.

The face value can be expressed as

The next upstream variable value is

And then the face value can be expressed as a normal interpolation by cell center values related with the face. All related information is included in the weight coefficientfβ.

in which

For a multi-dimensional unstructured grid as shown in Fig.4,Uα is not available, we can calculate Uα?as follows.

In order to ensure thatUα?is bounded, we assume that

From Eq.(20), it can be seen that the variable value at the interface is obtained from the central interpolation values of upstream and downstream.

This concept is firstly proposed by Onno Ubbink in 1997 in his doctoral thesis article and since then has received much attention[19]. An algorithm named Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAMS) is constructed in his thesis. This method is implemented in the in-house solver to simulate the sloshing in a swaying rectangular tank.

4. Simulations and discussions

4.1 Simulation model

The experimental data published by South Korea’s Daewoo Shipping and Marine Engineering Company in 2005 and used by Many authors[20], recommended by the 23rd ITTC Committee as a benchmark case[21], are used to verify the present method.

The rectangular tank used in the experiment is 0.8 m long, 0.5 m high and 0.35 m wide. The simulation is carried out with 30% filling height liquid of the tank. The dimensionless excitation frequency defined as the ratio between the excitation frequency and the liquid natural frequency, is 0.8, 0.9, 1.0, 1.2, respectively. The arrangement of the pressure monitoring points is shown in Fig.5, and sway displacement amplitudes are all 0.02 m. The sample frequency is 20 kHz in the experiments.

4.2 Analytical results

For a small amplitude wave, a linear analytical solution can be derived and used to verify the numerical results[22,23]. The wave height of the free surface can be formulated as follows[24]

Here a case of 30% filling height of liquid with excitation frequency ratio 0.8 and amplitude of 5×10-4m is simulated and the results are shown in Fig.6. The numerical results agree well with the analytical solution, which shows that the present method is reasonable and can be used to simulate the sway tank sloshing.

4.3 Comparison of 2-D and 3-D results

Firstly a 3-D real size model is used to simulate sloshing. The pressure time history is shown in Fig.7. A 3-D simulation is carried out with the same size as the experiment and a 2-D simulation is carried out of the middle plane in the span-wise direction. Pressures on four monitor points which cover areas both under and above the static water surface are obtained. The results show that the 2-D simulation is well consistent with the 3-D simulation and the 3-D effect on the pressure is negligible. Therefore, only 2-D simulation is carried out hereafter.

4.4 Grid independency

With the 2-D numerical model, the grid is generated with the pre-processor GAMBIT, a commercial software. Non-uniform quadrilateral grids with three different grid sizes 80×50, 160×100 and 240×150 are used to test the grid-independence. The grid is refined in a bi-exponential law with the ratio of 0.8. Time step is adjusted dynamically by a given courant number. As can be seen from Table 1, different grids do not affect the time history of the pressure, and affect its amplitude only slightly. However, the amount of computation will increase greatly with the grid refining. So 80×50 grids are used in this article.

4.5 Calculation results for 30% of filling height of liquid

The grid of 80×50 with 4 000 quadrilateral elements in total is used in the following simulations. Four different excitation frequencies are consideredfor 30% of filling height of liquid , which are 0.8, 0.9, 1.0, and 1.2. The results are shown in the following figures, in which, Exp is Daewoo’s experimental data, CFX is Ansys’ CFX11 results, and Gtea is this program’s calculation results.

It can be seen from Fig.8, for low-load conditions, the amplitude results of this method agree well with CFX’s experimental results. When the excitation frequency is low, the phase is in good agreement with the experiment. But, there is still a gap deviation near the low-level sloshing natural frequency.

4.6 Comparison with different filling height of liquid

The results with different filling heights of liquid with the same excitation frequency ratio ω/ωn=1.0 are shown in Fig.9. A high filling level corresponds a short period. This can be explained in this way. The high filling level means a high natural frequency or a high excitation frequency. From the pressure history at the monitor Points c1to c6, we can see that the Point c5is a transition point. Before it the pressure is low and wide, but the other points, on the other hand, see high and sharp impact.

4.7 Free surface shape

The snaps of the free surface in cases of different filling heights of liquid and the same excitation frequency are shown in Figs.10 through 12. Both CFX and GTEA results agree with the experiment. For the low filling height liquid case, the free surface is smooth, while in the high filling height liquid case, over-turning occurs. The present method can simulate this process well.

The time history of the wave height at three positions is shown in Fig.13, where it is seen that the numerical results are quite different from the analytical results which is anti-symmetrical with position and divergent[22,23]. This means that thenonlinear effect is important when the excitation frequency is near the resonant frequency. The natural frequency increases with the static liquid level and the period decreases. Numerical results are reasonable.

5. Conclusions

This article proposes a finite volume high order numerical algorithm based on NVD to simulate the sloshing on hybrid unstructured grids. The basic concept of this algorithm is that the free surface can be treated as a physical discontinuity, which can be captured by a well designed high order scheme. Non-inertial coordinate system is used for the rectangular box motion in order to avoid the moving mesh procedure. With the above mentioned methods, the sloshing in the sway tank can be simulated on unstructured fixed grids. This model is implemented in an in-house solver GTEA.

The experimental data published by South Korea’s Daewoo Shipping and Marine Engineering Company in 2005, recommended by the 23rd ITTC Committee as a benchmark case, are used to verify the present method. From the comparison of the pressure history on the monitor points, numerical results of both present GTEA and commercial software CFX agree well with the experiment data while GTEA shows some advantages on the secondary peak simulation. It is found that there is a transient zone for all simulation cases on one side (such as Point c1to c5) of which the impact is in a low and wide area, but the other side (such as Point c6) sees a high and sharp impact. The present method can simulate the free surface shape well and the wave height change is reasonable.

[1] PANIGRAHY P. K., SAHA U. K. and MAITY D. Experimental studies on sloshing behavior due to horizontal movement of liquids in baffled tanks[J]. Ocean Engineering, 2009, 36(3): 213-222.

[2] YAMAMOTO S., KATAOKA F. and SHIODA S. et al. Study on impact pressure due to sloshing in midsized LNG carrier[J]. International Journal of Offshore and Polar Engineering, 1995, 1(5): 10-16.

[3] TABRI K., MATUSIAK J. and VARSTA P. Sloshing interaction in ship collisions-An experimental and numerical study[J]. Ocean Engineering, 2009, 36(17-18): 1366-1376.

[4] FALTINSEN O. M., ROGNEBAKKE O. F. and LUKOVSKY I. A. et al. Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth[J]. Journal of Fluid Mechanics, 2000, 407: 201-234.

[5] HILL D. F. Transient and steady-state amplitudes of forced waves in rectangular basins[J]. Physics of Fluids, 2003, 15(6): 1576-1587.

[6] CELEBI M. S., KIM M. H. and BECK R. F. Fully non-linear 3-D numerical wave tank simulation[J]. Journal of Ship Research, 1998, 42(1): 33-45.

[7] KIM M. H., CELEBI M. S. and KIM D. J. Fully non-linear interactions of waves with a three-dimensional body in uniform currents[J]. Applied Ocean Research, 1998, 20(5): 309-321.

[8] LEE D. Y., CHOI H. S. Study on sloshing in cargo tanks including hydro elastic effects[J]. Journal of Marine Science and Technology, 1999, 4(1): 27-34.

[9] MORTON K., BAINES M. Numerical methods for fluid dynamics[M]. New York: Academic Press, 1982, 273-285.

[10] YANG Chi, L?HNER Rainald, LU Hai-dong. An unstructured grid based volume of fluid method for extreme wave and freely floating structure interactions[C]. Proceedings Conference of Global Chinese Scholars on Hydrodynamics. Shanghai, 2006, 415-422.

[11] AI Cong-fang, JIN Sheng. Three dimensional non-hydrostatic models for free surface flows with unstructured grid[J]. Journal of Hydrodynamics, 2008, 20(1): 108-116.

[12] YANG X. F., JAMES A. J. Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids[J]. Journal of computational physics, 2006, 214(1): 41-54.

[13] ZHENG Kun, SUN Zhao-chen and SUN Jia-wen. Numerical simulations of water wave dynamics based on SPH methods[J]. Journal of Hydrodynamics, 2009, 21(6): 843-850.

[14] REBOUILLAT S., LIKSONOV D. Fluid-structure interaction in partially filled liquid containers: A comparative review of numerical approaches[J]. Computers and Fluids, 2010, 39(5): 739-746.

[15] IBRAHIM R., PILIPCHUK V. and IKEDA T. Recent advances in liquid sloshing dynamics[J]. Applied Mechanics Review, 2001, 54(2): 133-199.

[16] LEONARD B. P. The ULTIMATE conservative difference scheme applied to unsteady one dimensional advection[J]. Computational Method in Applied Mechanics and Engineering, 1991, 88(1): 17-74.

[17] DARWISH M. S., MOUKALLED F. TVD schemes for unstructured grids[J]. International Journal of Heat and Mass Transfer, 2003, 46(4): 599-611.

[18] Ubbink Onno. Numerical prediction of two fluid systems with sharp interfaces[D]. Ph. D. Thesis, London, UK: University of London 1997.

[19] TAO Wen-quan. Modern progress on computational heat transfer[M]. Beijing: Science Press, 2001, 105-158(in Chinese).

[20] CHEN Y. G., PRICE W. G. Numerical simulation of liquid sloshing in a partially filled container with inclusion of compressibility effects[J]. Physics of Fluids, 2009, 21(11): 112105.

[21] THE ITTC Seakeeping Committee. Final report and recommendation to the 23rd ITTC[C]. The 23rd proceedings of International Towing Tank conference. Venice, Italy, 2002,1-5.

[22] FALTINSEN O. M. A numerical non-linear method for sloshing in tanks with two dimensional flow[J]. Journal of Ship Research, 1978, 18(4): 224-241.

[23] WU G. X. Second order resonance of sloshing in a tank[J]. Ocean Engineering, 2007, 34(17-18): 2345-2349.

[24] ZHENG Xing. Analysis and improvement of algorithm for SPH and application for free surface flow simulation[D]. Harbin: Harbin Engineering University, 2009(in Chinese).

10.1016/S1001-6058(09)60126-8

* Project supported by the China Postdoctoral Science Foundation (Grant No. 20100471016), the Fundamental Research Funds for Major Universities (Grant No. HEUCF 100307).

Biography: MING Ping-jian (1980-), Male, Ph. D., Lecturer