Numerical Simulation of Nonlinear Sloshing Waves in Three-dimensional Tank based on DBIEM

XÜ Gang,MA Xiao-jian,LIU Yong-tao,ZHU Ren-qing

(School of Naval Architecture and Ocean Engineering,Jiangsu University of Science and Technology,Zhenjiang 212003,China)

Numerical Simulation of Nonlinear Sloshing Waves in Three-dimensional Tank based on DBIEM

XÜ Gang,MA Xiao-jian,LIU Yong-tao,ZHU Ren-qing

(School of Naval Architecture and Ocean Engineering,Jiangsu University of Science and Technology,Zhenjiang 212003,China)

Based on the fully nonlinear velocity potential theory,the liquid sloshing in a three-dimensional tank is studied.The governing Laplace equation with fully nonlinear boundary conditions on the moving free surface is solved by using the indirect desingularized boundary integral equation method(DBIEM).The fourth-order predictor-corrector Adams-Bashforth-Moulton scheme(ABM4) and mixed Eulerian-Lagrangian(MEL)method are used for the time-stepping integration of the free surface boundary conditions.A smoothing scheme,B-spline curve,is applied to both the longitudinal and transverse directions of the tank to eliminate the possible saw-tooth instabilities.When the tank is undergoing horizontal regular motion of small amplitude,the calculated results are in very good agreement with linear analytical solution.

sloshing waves;DBIEM;fully nonlinear

0 Introduction

Sloshing motions of liquid are associated with various engineering problems,such as liquid cargo in a LNG carrier,the liquid oscillations in large storage tanks caused by earthquakes, the motions of liquid fuel in aircraft and spacecraft,the liquid motions in containers and the water flow on the decks of ships.As the motion can become large or even violent when resonance occurs,the liquid load can cause structural damage on the container and/or lead to loss of stability of the liquid carrier such as a ship.

Since sloshing has important implications in engineering,it has been extensively studied over the past years.Abramson[1]used a linear theory to simulate small amplitude sloshing in a container,Solaas and Faltinsen[2]employed a perturbation theory.Wu[3]considered the second-order resonance conditions in a rectangular tank.Firouz-Abadi et al[4]used BEM model to investigate second-order analysis of sloshing in tanks with arbitrary shapes under both horizontal and vertical excitations.For large amplitude sloshing,Faltinsen[5]adopted the boundaryelement method to simulate fluid motion in a 2D rectangular container subjected to a horizontal excitation.Cho and Lee[6]presented a non-linear finite element method for the simulation of large amplitude sloshing in a rectangular baffled tank,subjected to horizontal forced excitation.They performed parameter studies on the effects of the baffle on the non-linear liquid sloshing.Chen and Nokes[7]applied Navier-Stokes equations to study 2D sloshing motion (surge,heave and pitch)in a tank via coordinate transformation.Wu et al[8]used the finite element method(FEM)and analyzed 3D sloshing waves through the fully nonlinear velocity potential theory.Liu and Lin[9]used the spatially averaged Navier-Stokes equations to study 3D nonlinear liquid sloshing.Frandsen[10]conducted a series of numerical experiment in a 2D tank which is moved both horizontally and vertically via σ-coordinate transformation.Wu and Chen[11]used 3D finite difference method to solve the wave sloshing in a 3D tank excited by coupled surge and sway motions.For random sloshing,Wang and Khoo[12]adopted FEM and considered the 2D nonlinear sloshing problem in a tank under random excitation.Sriram et al[13]performed the finite element analysis of non-linear sloshing in a rectangular tank under both horizontal and vertical random excitations.

In the present study,the desingularized boundary integral equation method(DBIEM)is selected for simulating 3D sloshing wave.The method has been successfully used previously in solving nonlinear water wave problems,such as in the work by Beck[14],Kim et al[15],Celebi[16], Kara et al[17]and Zhang et al[18].The main advantage of DBIEM,compared with FEM,lies in having only to discretize the surface of fluid domain.When the boundary of the fluid domain is confined and the number of the discretized elements is limited,the DBIEM may offer a better computational efficiency and less memory requirement,even its matrix is fully populated.Compared with the conventional BEM,The integral kernels of the DBIEM are no longer singular as the singularities are placed slightly outside the fluid domain.This is particularly advantageous when the direct differentiation is applied to the integral equation to obtain the velocity.

In this study,we shall focus on the 3D sloshing problem in a rectangular tank.Similar problems have been considered by Wu et al[8]and Kara et al[17].Here the DBIEM is employed to solve the boundary value problem at each time step.The fourth-order predictor-corrector ABM4 scheme and mixed Eulerian-Lagrangian(MEL)method are used for the time-stepping integration of the free surface boundary conditions.Since wave breaking is not considered in this work,the position of the nodes on free surface is tracked by applying semi-Lagrangian approach(Zhang et al[19];Zhang et al[20];Khoo and Kim[21]),in which the nodes on free surface are allowed to move only in vertical direction,with the horizontal motion of the nodes on the free surface held fixed.This approach has the advantage of avoiding the task of re-gridding the free surface at each time step.For stable time-step simulation,a B-spline smoothing scheme is applied in both longitudinal and transverse directions of the tank to prevent saw-tooth instability,and the smoothing scheme is used at every three time steps.Numerical results obtained show that the present desingularized model is effective in the simulation for 3D sloshing waves.

1 Mathematical formulation

A Cartesian co-ordinate system oxyz is defined for 3D sloshing waves.The origin is at the centre of the undisturbed free surface,as shown in Fig.1.x and y are in the longitudinal and transverse directions of the rectangular tank,z points vertically upwards.

Based on potential flow theory,the velocity potential φ in computational fluid domain D satisfies the Laplace equation:

Fig.1 The co-ordinate system and the sloshing tank

where U is the velocity of the tank and n is outward normal vector of the surfaces of the tank wall and bottom.

On the instantaneous free surface ΓF,the dynamic and kinematic conditions can be written as:

where t and g denote time and gravitational acceleration,respectively.The initial conditions may be expressed as

2 Desingularized boundary integral equation method

In this study,the indirect DBIEM is employed to solve the boundary value problem for the unknown velocity potential φ（ x,y,z,t）at each time step.This method obtains the solution by distributing Rankine sources over a surface S outside the fluid domain D.This surface is at a small distance away from the corresponding real boundary of the fluid.The velocity potential in the fluid domain D can be written as follows:

For the problem considered in this work,we construct the solution using a constantstrength source point within each element over the integration boundary SFand a constantstrength source point over the integration surface SW,where SFis the integration surface above the free surface ΓF,and SWis the integration surface outside the real boundary ΓWof the tank. That is

By applying the boundary conditions,we obtain boundary integral equations for the unknown strength of the singularities,σF（q,t）and σW（q,t）,respectively:

In the desingularized method,the source distribution is outside the fluid domain so that the source points never coincide with the field points(control or collocation points)and therefore the integrals are non-singular.In addition,because of the desingularization,we can use simple isolated Rankine sources and obtain the equivalent accuracy.This greatly reduces the complexity of the form of the influence coefficients that make up the elements of the kernel matrix(Zhang et al[19]).Then the integral equations in Eq.(10)and Eq.(11)can be replaced by a discrete summation of N-isolated singularities located at a small distance away from the corresponding nodal point on the free surface and body,

The desingularized distance between isolated source point and corresponding nodal point is given by

where ldand β are constants and Dmis a measure of the local mesh size(typically the square root of the local mesh area).The accuracy and convergence of the solutions are sensitive to the choices of ldand β.Therefore,appropriate ldand β values need to be determined after numerical test.The recommended values are ld=0.5-1.0 and β=0.5.A detailed study with regard tothe performance of DBIEM with the desingularization parameters was reported by Cao et al22].

Once the above integral equations using isolated Rankine source are solved at each time step,the fluid velocity in Eq.(3)and Eq.(4)can be calculated from direct derivatives,

3 Time-stepping integration scheme

In order to obtain the velocity potential and free surface elevation at each time step,the fourth-order predictor-corrector Adams-Bashforth-Moulton scheme(ABM4)and mixed Eulerian-Lagrangian(MEL)method are used.Integrating Eq.(3)and Eq.(4)by ABM4 and MEL is called time marching.Using the total derivative δ/δt=∂/∂t+v→·▽,the fully nonlinear free surface conditions can be modified as follows in Lagrangian frame,

In the semi-Lagrangian approach,a time-stepping integration procedure must be employed to obtain the values of velocity potential and wave elevation on the instantaneous free surface.After solving the boundary value problem and obtaining the fluid velocity on the free surface at each time step,the free surface boundary conditions can be treated as ordinary differential equations to be marched in time.The general form of the dynamic and kinematic boundary conditions Eq.(18)and Eq.(19)can be rewritten as

ABM4 scheme(Zhang et al[19])is selected for integrating Eq.(20)and Eq.(21)with time.It is a fourth-order predictor and corrector method which only requires two evaluations of the functions g（ η,φ,t）and f（ η,φ,t）at each time step.In the ABM4 scheme,the velocity potential and wave elevation are firstly predicted by Adams-Bashforth method as follows:

and then these are iteratively corrected by Adams-Moulton algorithm,

where△t is the time step.

4 Linearised solution of wave elevation

For the periodic oscillation,i.e.U=Aωsinωt,Wu[3]gave the linearized analytical solution for wave elevation η:

where kn=nπ/L,L is the length of the tank.If only the motion of the tank in x direction will be considered,then

If we further consider the problem with coupled surge and sway motions,i.e.Ux=Axωxsinωxt and Uy=Ayωysinωyt,the 3D linearized analytical solution of free surface elevation can be given as follows:

where B is the breadth of the tank,

It should be noted that the natural frequencies in the 3D cases are(Liu and Lin,2008):

5 Numerical results and discussions

5.1 One-directional periodic oscillation

We first consider one directional periodic oscillation with

Fig.2 Mesh on the surface of rectangular tank(lower),mesh on the free surface(upper)

where ω is excitation frequency.

In the study,the dimensions of the tank are chosen as L=1.0 m,B=0.5 m and h=0.5 m where L,B and h are the length,the width and the water depth,respectively.A typical initial mesh is illustrated in Fig.2.

Fig.3 Comparison of free surface elevation at(-0.476 2 m,0,0)in an excited tank between the present fully nonlinear solution(red dotted line)and analytical solution(black line)

Fig.4 Free surface profiles for ωx=0.999ω0x(t=2.5～25.0 s,interval is 5 s) due to one-directional motion

5.2 Two-directional motions

The case with the motions in both x direction and y direction has been considered.ωx= 0.5ω0x,ωy=0.5ω0yand ωx=0.9ω0x,ωy=0.9ω0yare chosen as excitation frequency,where ω0xand ω0yare the lowest frequency of fluid in the tank due to the x-directional and y-directional motion,respectively.Ax=Ay=0.01 m and Ax=Ay=0.000 5 m are chosen as amplitude in surge and sway modes,respectively.Fig.5 presents the time history of free surface elevation at (-0.476 2 m,-0.227 3 m,0).Comparison between analytical solution and numerical results shows that they are also in good agreement.Fig.6 shows the snapshots of the free surface pro-files between t=2.5 s and t=25.0 s,at the intervals equal to 2.5 s,for the case of ωx=0.9ω0x, ωy=0.9ω0yand Ax=Ay=0.000 5 m.

Fig.5 Comparison of free surface elevation at(-0.476 2 m,-0.227 3 m,0)in an excited tank between the present fully nonlinear solution(red dotted line)and analytical solution(black line)

Fig.6 Free surface profiles at different time(t=2.5～25.0 s,interval is 5 s)due to two-directional motions(ωx=0.9ω0x,ωy=0.9ω0y,Ax=Ay=0.000 5 m)

5.3 Three-directional motions

The case with the motions in x,y and z directions has been considered finally.ωx=0.5ω0x, ωy=0.5ω0y,ωz=0.5ω0xand ωx=0.9ω0x,ωy=0.9ω0y,ωz=0.9ω0xare chosen as excitation frequency,where Ax=Ay=Az=0.01 m and Ax=Ay=Az=0.000 5 m are chosen as amplitude in surge, sway and heave modes,respectively.Fig.7 presents the time history of free surface elevation at the four corners of the tank.

Fig.7 Wave elevation history at four corners

6 Conclusions

In this paper,DBIEM coupled with MEL time marching scheme is applied to simulate sloshing waves in a 3D tank undergoing specified horizontal motion.The fourth-order predictor-corrector ABM4 scheme is used for the time-stepping integration of the free surface bound-ary conditions.The position of instantaneous free surface is tracked by applying semi-Lagrangian approach.The saw-tooth instability is overcome by applying B-spline smoothing scheme to both the longitudinal and transverse directions during the simulation.The model is validated by available linear theory and checking the conversation of fluid mass.All the numerical results agree fairly with the linear analytical solution for small amplitude cases.It is our interest in future works that the present model is extended to development of 3D sloshing wave in Liquefied Natural Gas Carrier.

[1]Abramson H N.The dynamic behavior of liquid in moving containers[R].NASA Report,SP 106,1996.

[2]Solaas F,Faltinsen O M.Combined numerical solution for sloshing in two-dimensional tanks of general shape[J].J Ship Res.,1997,41:118-129.

[3]Wu G X.Second-order resonance of sloshing in a tank[J].Ocean Eng.,2007,34:2345-2349.

[4]Firouz-Abadi R D,Ghasemi M,Haddadpour H.A modal approach to second-order analysis of sloshing using boundary element method[J].Ocean Eng.,2011,38:11-21.

[5]Faltinsen O M.A numerical non-linear method of sloshing in tanks with two-dimensional flow[J].J Ship Res.,1978,22: 193-202.

[6]Cho J R,Lee H W.Numerical study on liquid sloshing in baffled tank by nonlinear finite element method[J].Comput. Methods.Appl.Mech.Eng.,2004,193:2581-2598.

[7]Chen Bang-Fuh,Nokes R.Time-independent finite difference analysis of fully non-linear and viscous fluid sloshing in a rectangular tank[J].J Comput.Phys.,2005,209:47-81.

[8]Wu G X,Ma Q W,Eatock Taylor R.Numerical simulation of sloshing waves in a 3D tank based on a finite element method[J].Appl.Ocean Res.,1998,20:337-355.

[9]Liu D M,Lin P Z.A numerical study of three-dimensional liquid sloshing in tanks[J].J Comput.Phys.,2008,227: 3921-3939.

[10]Frandsen Jannette B.Sloshing motions in excited tanks[J].J Comput.Phys.,2004,196:53-87.

[11]Wu Chih-Hua,Chen Bang-Fuh.Sloshing waves and resonance modes of fluid in a 3D tank by a time-independent finite difference method[J].Ocean Eng.,2009,36:500-510.

[12]Wang C Z,Khoo B C.Finite element analysis of two-dimensional nonlinear sloshing problems in random excitations[J]. Ocean Eng.,2005,32:107-133.

[13]Sriram V,Sannasiraj S A,Sundar V.Numerical simulation of 2D sloshing waves due to horizontal and vertical random excitation[J].Appl.Ocean Res.,2006,28:19-32.

[14]Beck R F.Time-domain computations for floating bodies[J].Appl.Ocean Res.,1994,16:267-282.

[15]Kim M H,Celebi M S,Kim D J.Fully nonlinear interactions of waves with a three-dimensional body in uniform currents[J].Appl.Ocean Res.,1998,20:309-321.

[16]Celebi M S.Nonlinear transient wave-body interactions in steady uniform currents[J].Comput.Methods.Appl.Mech. Eng.,2001,190:5149-5172.

[17]Kara F,Tang C H,Vassalors D.Time domain three-dimensional fully nonlinear computations of steady body-wave interaction problem[J].Ocean Eng.,2007,34:776-789.

[18]Zhang X S,Bandyk P,Beck Robert F.Seakeeping computations using double-body basis flows[J].Appl.Ocean Res., 2010,32:471-482.

[19]Zhang X T,Khoo B C,Lou J.Wave propagation in a fully nonlinear numerical wave tank:a desingularized method[J]. Ocean Eng.,2006,33:2310-2331.

[20]Zhang X T,Khoo B C,Lou J.Application of desingularized approach to water wave propagation over three-dimensional topography[J].Ocean Eng.,2007,34:1449-1458.

[21]Koo W C,Kim M H.Fully nonlinear wave-body interactions with surface-piercing bodies[J].Ocean Eng.,2007,34: 1000-1012.

[22]Cao Y,Schultz W W,Beck R F.Three dimensional desingularized boundary integral methods for potential problems[J]. Int.J Numer.Methods.Fluids,1991,12:785-803.

基于无奇异边界元法模拟三维全非线性液舱晃荡

徐刚，马小剑，刘永涛，朱仁庆
（江苏科技大学船舶与海洋工程学院，江苏镇江212003）

文章基于全非线性势流理论对三维液舱晃荡进行了数值模拟，其控制方程由无奇异边界积分方程法（Desingularized Boundary Integral Equation Method,DBIEM）进行离散求解，在求解全非线性的自由面微分方程时，文中采用混合欧拉—拉格朗日法（Mixed Eulerian-Lagrangian,MEL）和四阶Adams-Bashforth-Moulton（ABM4）预报—修正方法，为了避免结果发散即增强数值稳定性，文中采用B样条法来光顺自由面。在微幅水平激励下，该文中得到的结果与解析解吻合较好。

声辐射模态；单层阵列；声场分离技术；近场声全息

O35

:A

国家自然科学基金资助（51309125,51409128,51379094,51179077）；

徐刚（1981-），男，博士，江苏科技大学船舶与海洋工程学院副教授；

O35

:A

10.3969/j.issn.1007-7294.2017.06.002

1007-7294（2017）06-0661-11

江苏高校优势学科建设工程资助项目资助

马小剑（1982-），男，博士，江苏科技大学船舶与海洋工程学院讲师；

date：2016-12-28

Supported by the National Natural Science Foundation of China(Grant Nos.51309125, 51409128,51379094,51179077)and the Project Founded by Priority Academic Program Development of Jiangsu Higher Education Institutions

Biography：XU Gang(1981-),male,Ph.D.associate prof.,E-mail:me_xug@qq.com;

MA Xiao-jian(1982-),make,Ph.D.,lecturer.

刘永涛（1977-），男，博士，江苏科技大学船舶与海洋工程学院副教授；

朱仁庆（1965-），男，博士，江苏科技大学船舶与海洋工程学院教授。