# Nonlinear Finite Element Analysis of Sloshing

Authors: Siva Srinivas Kolukula; P. Chellapandi
Journal: Advances in Numerical Analysis
 (2013)

---

## Abstract

The disturbance on the free surface of the liquid when the liquid-filled tanks are excited is called sloshing. This paper examines the nonlinear sloshing response of the liquid free surface in partially filled two-dimensional rectangular tanks using finite element method. The liquid is assumed to be inviscid, irrotational, and incompressible; fully nonlinear potential wave theory is considered and mixed Eulerian-Lagrangian scheme is adopted. The velocities are obtained from potential using least square method for accurate evaluation. The fourth-order Runge-Kutta method is employed to advance the solution in time. A regridding technique based on cubic spline is employed to avoid numerical instabilities. Regular harmonic excitations and random excitations are used as the external disturbance to the container. The results obtained are compared with published results to validate the numerical method developed.


---

## Body

## 1. Introduction

It is common everyday knowledge to each of us that any small container filled with liquid must be moved or carried very carefully to avoid spills. For example, one has to be careful while carrying a cup of coffee while moving, because the motion of the person makes coffee spill. Such a motion on the free surface of the liquid, due to external excitation in the liquid-filled containers, is called sloshing. Sloshing is likely to be seen whenever we have a liquid with a free surface in the presence of gravity. At equilibrium the free surface of the liquid is static and coincides with a gravitational equipotential surface. When the surface is perturbed, an oscillation is set up in which the energy oscillates between kinetic energy and gravitational potential energy. The phenomenon called sloshing occurs in a variety of engineering applications such as sloshing in liquid-propellant launch vehicles, sloshing in liquids used in industries to store oil, water, chemicals, liquefied natural gases, and so forth, and sloshing in the nuclear reactors of pool type, nuclear fuel storage tanks under earthquake. The liquid sloshing may cause huge loss of human, economic, and environmental resources owing to unexpected failure of the container; for example the spillage of toxic chemicals stored in tanks in industries can cause contamination of soil and the environment. Thus, understanding the dynamic behaviour of liquid free surface is essential. As a result, the problem of sloshing has attracted many researchers and engineers targeting to understand the complex behaviour of sloshing and to design the structures to withstand its effects.Abundant research has been made on the sloshing phenomenon and the literature is vast with wide varieties of numerical methods, analytical solutions, and experiments. The position of the free surface of liquid is not known a priori and obeys a dynamic boundary condition which is nonlinear and makes the problem of sloshing that a nonlinear boundary value problem. Although the sloshing problem is nonlinear, by assuming the free-surface elevation to be small and applying the linearized free-surface boundary condition, a linear theory of sloshing is developed [1, 2]. This linear theory is acceptable in few cases when the amplitude of external excitation is small and not in the neighborhood of the sloshing natural frequency. When the above mentioned conditions do not hold, linear theory fails to predict the sloshing response accurately. Hence nonlinear analysis becomes inevitable for accurate and reliable evaluation of liquid sloshing.Nonlinear sloshing problem is difficult to solve analytically because of its nonlinear boundary conditions implying a numerical modeling is necessary. Faltinsen [3, 4] solved the nonlinear sloshing problem numerically and derived the analytical solution using perturbation approach with two-dimensional flow. Nakayama and Washizu [5] employed the boundary element method for the problem. Wu and Taylor [6, 7] applied finite element analysis for two dimensional nonlinear transient waves. Chen et al. [8] applied finite difference method to simulate large-amplitude sloshing under seismic excitations. Turnbull et al. [9] used sigma-transformed finite element inviscid flow solver for the problem. Frandsen [10] analysed the nonlinear sloshing motions of liquid under vertical, horizontal, and combined motions of the tank using analytical and numerical methods. The author used perturbation methods for analytical solution and modified sigma-transformed finite difference method for numerical solution.  Cho and Lee [11] used semi-Lagrangian nonlinear finite element approximation to analyse the large amplitude sloshing in two-dimensional tanks. Wang and Khoo [12] analysed the nonlinear sloshing in two-dimensional tanks under random excitations. Sriram et al. [13] analysed nonlinear sloshing in two-dimensional tanks using finite element and finite difference method.  Biswal et al. [14] analysed the nonlinear sloshing response in tanks with baffles using finite element method. Ibrahim et al. [15] give an excellent review of sloshing phenomenon with extensive number of references available in the literature.In the present paper a numerical approach based on mixed Eulerian-Lagrangian scheme is adopted. The free surface nodes behave like Lagrangian particles and interior nodes behave like Eulerian particles. The nonlinear sloshing analysis is carried out using finite element method. A four-noded isoparametric element is used in the analysis. The calculation of velocities from velocity potential is an important step to study the sloshing behaviour. Thus velocities must be calculated accurately for accurate sloshing analysis. The velocity field is interpolated from the velocity potential according to least square method [16]. Fourth-order Runge-Kutta method is employed to advance the solution in time. As the time proceeds in the simulation due to Lagrangian behaviour of the free surface nodes, these nodes move closer and develop a steep gradient leading to numerical instability. To get rid of this problem, a cubic spline interpolation is used for the regridding the free surface uniformly. In the present simulation the tank is assumed to be rigid with aspect ratio (h/L) of 0.5; h is depth of fluid, L is length of the tank. Sloshing response is simulated when the external excitation frequency is in resonance and off resonance region.  For horizontal excitations the free surface undergoes resonance when excitation frequency is equal to fundamental slosh frequency and shows beating phenomenon when excitation frequency is close to fundamental slosh frequency. The present numerical model is validated with Frandsen [10] numerical results for harmonic excitations and then applied for sloshing response due to random excitation.


## 2. Governing Equations

Consider a rectangular tank fixed in Cartesian coordinate systemOxz, which is moving with respect to inertial Cartesian coordinate system O0x0z0. The origins of this system are at the left end of the tank wall at the free surface and pointing upwards in z direction. These two Cartesian systems coincide when the tank is at rest. Figure 1 shows the tank in the moving Cartesian coordinate system Oxz along with the prescribed boundary conditions. The tank is assumed to get displaced along x-axes, and the tank position at time is given by
(1)X=Xt(t).
Fluid is assumed to be inviscid, incompressible, and irrotational. Therefore the fluid motion is governed by Laplace’s equation with the unknown as velocity potential ϕ:
(2)∇2ϕ=0.
The fluid obeys Neumann boundary conditions at the walls of the container and Dirichlet boundary condition at the liquid free surface. In the moving coordinate system the velocity component of the fluid normal to the walls is zero. Hence, on the bottom and on the walls of the tank (ΓB) we have
(3)∂ϕ∂n|x=0,L=0,∂ϕ∂n|z=-h=0.
On the free surface (Γs) dynamic and kinematic conditions hold; they are given as
(4)∂ϕ∂t|z=ζ+12∇ϕ·∇ϕ+gζ+xXt′′,∂ζ∂t+∂ϕ∂x∂ζ∂x-∂ϕ∂z=0.
Rewrite the above equations in the Lagrangian form [16]:
(5)dϕdt|z=ζ=12∇ϕ·∇ϕ-gζ-xXt′′,(6)dxdt=∂ϕ∂x,dzdt=∂ϕ∂z,
where ζ is the free surface elevation measured vertically above still water level, Xt′′ is the horizontal acceleration of the tank, and g is the acceleration due to gravity.Figure 1
Sloshing wave tank in moving coordinate system.Equations (1)–(6) give the complete behaviour of nonlinear sloshing in fluids. The position of the fluid free surface is not known a priori; to solve the problem the fluid is assumed to be at rest. Thus the initial conditions for the free surface in the moving Cartesian system at t=0 and z=0 can be written as:
(7)ϕ(x,0,0)=-xdXt(t)dt,(8)ζ(x,0)=0.
Using these initial conditions, Laplace equation (2) is solved and the free surface elevation and potential are updated for the subsequent time steps using (5)-(6).


## 3. Numerical Procedure

### 3.1. Finite Element Formulation

The solution of the nonlinear sloshing boundary value problem is obtained using finite element method. The entire liquid domain is discretized by using four-noded isoparametric quadrilateral elements. A typical mesh structure is shown in Figure2. By introducing the finite element shape functions the liquid velocity potential can be approximated as
(9)ϕ(x,z)=∑j=1nNj(x,z)ϕj,
where Nj is the shape function, n is the number of nodes in the element, and ϕj is nodal velocity potential.Figure 2
A typical mesh of the liquid domain using isoparametric four-noded elements.On applying Galerkin residual method to the Laplace equation, we get(10)Kϕ=0,
with matrix K defined by
(11)K=∫Ω(∇N)T(∇N)dΩ.
Matrix K is analogous to stiffness matrix in structural problems. Using (10), from the known free-surface velocity potential, the velocity potential in the interior nodes is calculated.


### 3.2. Velocity Recovery

To track the free surface (4) need velocities, which are computed from the calculated potential field using
(12)ν-=∇·ϕj.
The velocities calculated using (12) are the velocities at the Gauss integration points and they do not possess interelement continuity and have a low accuracy at nodes and element boundaries. Utmost care should be taken to calculate the velocities; a small error in the velocity recovery will affect the accuracy of free surface updating and gets accumulated with time and leads to underestimation of sloshing response. In order to derive a smoothed and continuous velocity, patch recovery technique [17] is applied. In patch recovery technique, the continuous velocity field is obtained by considering the linear interpolation of the velocities at the Gauss integration points:
(13)ν^=a1+a2ξ+a3η+a4ξη,
where ν^ is any velocity component (ν^x or ν^y), ξ, η are the Gauss locations, and a1,a2,a3,a4 are unknowns which need to be evaluated. To evaluate these unknowns, a least square fit is considered between ν^ and ν-:
(14)F(a)=∑i=14{ν^(ξi,ηi)-ν-(ξi,ηi)}2,
where i is 2×2 order Gauss integration points. Then, the four unknown coefficients are determined from four simultaneous equations obtained from
(15)∂F(a)∂ak=0,k=1,2,3,4.
Substituting the obtained ak’s in (13) gives the velocity values for individual elements and these are averaged for the common nodes. Finally, a smoothed velocity field which is interelement continuous is constructed by interpolating the finite element shape functions used in (9) and nodal averaged velocities. The global continuous velocity field is given as
(16)ν=N·ν^,
where ν is velocity component (νx or νy). Using the patch recovery technique velocity components νx and νy are calculated.


### 3.3. Numerical Time Integration and Free Surface Tracking

After calculating the velocity at a time stept, we need to calculate the position of free surface from (6) and determine the potential on the free surface using (5) for the next time step t+∇t. As a result, the liquid mesh and the boundary condition required for the next-time step are established. This is done using fourth-order Runge-Kutta explicit time integration method. The nodal coordinates of the free surface and the associated velocity potential at a current time step i are known and can be represented in a single variable as
(17)si=(xi,zi,ϕi),
where
(18)xi=(x1,x2,…,xNX+1)i,zi=(z1,z2,…,zNX+1)i,xi=(ϕ1,ϕ2,…,ϕNX+1)i,
where NX is number of segments along the free surface. Similarly the time derivative can be written as
(19)DsiDt=F(ti,si)=Fi.
The free surface position and associated velocity potential at the next time step i+1 can be expresses as
(20)si+1=si+s16+s23+s33+s46,
where
(3.3)s1=∇tF(ti,si),s2=∇tF(ti+∇t2,si+s12),s3=∇tF(ti+∇t2,si+s22),s4=∇tF(ti+∇t,si+s3).
After obtaining the new positions and potential of the free surface, the liquid domain is remeshed based on these obtained new coordinate positions.


### 3.4. Regridding Algorithm

At the beginning of the numerical simulation, the free surface nodes are uniformly distributed along thex-direction with zero surface elevation. As the time proceeds the free surface nodes are spaced unequally and cluster into a steep gradient leading to numerical instability. This problem occurs for a long time simulation; to avoid this instability an automatic regridding condition using cubic spline is employed when the movement of the nodes is 75% more or less then the initial grid spacing. For the regridding, first the free surface length Lf is obtained. Then the free surface is divided into NX segments with the identical arc length. The coordinates of node is denoted as (xl,yl)(l=1,2,…,NX+1) and let the arc length between two successive points l and l+1 be Sl. Being a uniform regridding, Sl can be expressed as
(22)Sl=lLfNX.
The coordinates of the nodes (xl,yl) is a function of the arc length Sl:
(23)(xl,yl)=f(Sl).
The cubic spline interpolation is used to calculate the coordinates (xl,yl) and the velocity potential on the new uniform free surface is also obtained in a similar fashion.


### 3.5. Complete Algorithm for Nonlinear Sloshing

Including all the steps above, the algorithm for numerical simulation of nonlinear sloshing is as shown in Figure3.Figure 3
Numerical procedure for nonlinear sloshing simulation.


## 3.1. Finite Element Formulation

The solution of the nonlinear sloshing boundary value problem is obtained using finite element method. The entire liquid domain is discretized by using four-noded isoparametric quadrilateral elements. A typical mesh structure is shown in Figure2. By introducing the finite element shape functions the liquid velocity potential can be approximated as
(9)ϕ(x,z)=∑j=1nNj(x,z)ϕj,
where Nj is the shape function, n is the number of nodes in the element, and ϕj is nodal velocity potential.Figure 2
A typical mesh of the liquid domain using isoparametric four-noded elements.On applying Galerkin residual method to the Laplace equation, we get(10)Kϕ=0,
with matrix K defined by
(11)K=∫Ω(∇N)T(∇N)dΩ.
Matrix K is analogous to stiffness matrix in structural problems. Using (10), from the known free-surface velocity potential, the velocity potential in the interior nodes is calculated.


## 3.2. Velocity Recovery

To track the free surface (4) need velocities, which are computed from the calculated potential field using
(12)ν-=∇·ϕj.
The velocities calculated using (12) are the velocities at the Gauss integration points and they do not possess interelement continuity and have a low accuracy at nodes and element boundaries. Utmost care should be taken to calculate the velocities; a small error in the velocity recovery will affect the accuracy of free surface updating and gets accumulated with time and leads to underestimation of sloshing response. In order to derive a smoothed and continuous velocity, patch recovery technique [17] is applied. In patch recovery technique, the continuous velocity field is obtained by considering the linear interpolation of the velocities at the Gauss integration points:
(13)ν^=a1+a2ξ+a3η+a4ξη,
where ν^ is any velocity component (ν^x or ν^y), ξ, η are the Gauss locations, and a1,a2,a3,a4 are unknowns which need to be evaluated. To evaluate these unknowns, a least square fit is considered between ν^ and ν-:
(14)F(a)=∑i=14{ν^(ξi,ηi)-ν-(ξi,ηi)}2,
where i is 2×2 order Gauss integration points. Then, the four unknown coefficients are determined from four simultaneous equations obtained from
(15)∂F(a)∂ak=0,k=1,2,3,4.
Substituting the obtained ak’s in (13) gives the velocity values for individual elements and these are averaged for the common nodes. Finally, a smoothed velocity field which is interelement continuous is constructed by interpolating the finite element shape functions used in (9) and nodal averaged velocities. The global continuous velocity field is given as
(16)ν=N·ν^,
where ν is velocity component (νx or νy). Using the patch recovery technique velocity components νx and νy are calculated.


## 3.3. Numerical Time Integration and Free Surface Tracking

After calculating the velocity at a time stept, we need to calculate the position of free surface from (6) and determine the potential on the free surface using (5) for the next time step t+∇t. As a result, the liquid mesh and the boundary condition required for the next-time step are established. This is done using fourth-order Runge-Kutta explicit time integration method. The nodal coordinates of the free surface and the associated velocity potential at a current time step i are known and can be represented in a single variable as
(17)si=(xi,zi,ϕi),
where
(18)xi=(x1,x2,…,xNX+1)i,zi=(z1,z2,…,zNX+1)i,xi=(ϕ1,ϕ2,…,ϕNX+1)i,
where NX is number of segments along the free surface. Similarly the time derivative can be written as
(19)DsiDt=F(ti,si)=Fi.
The free surface position and associated velocity potential at the next time step i+1 can be expresses as
(20)si+1=si+s16+s23+s33+s46,
where
(3.3)s1=∇tF(ti,si),s2=∇tF(ti+∇t2,si+s12),s3=∇tF(ti+∇t2,si+s22),s4=∇tF(ti+∇t,si+s3).
After obtaining the new positions and potential of the free surface, the liquid domain is remeshed based on these obtained new coordinate positions.


## 3.4. Regridding Algorithm

At the beginning of the numerical simulation, the free surface nodes are uniformly distributed along thex-direction with zero surface elevation. As the time proceeds the free surface nodes are spaced unequally and cluster into a steep gradient leading to numerical instability. This problem occurs for a long time simulation; to avoid this instability an automatic regridding condition using cubic spline is employed when the movement of the nodes is 75% more or less then the initial grid spacing. For the regridding, first the free surface length Lf is obtained. Then the free surface is divided into NX segments with the identical arc length. The coordinates of node is denoted as (xl,yl)(l=1,2,…,NX+1) and let the arc length between two successive points l and l+1 be Sl. Being a uniform regridding, Sl can be expressed as
(22)Sl=lLfNX.
The coordinates of the nodes (xl,yl) is a function of the arc length Sl:
(23)(xl,yl)=f(Sl).
The cubic spline interpolation is used to calculate the coordinates (xl,yl) and the velocity potential on the new uniform free surface is also obtained in a similar fashion.


## 3.5. Complete Algorithm for Nonlinear Sloshing

Including all the steps above, the algorithm for numerical simulation of nonlinear sloshing is as shown in Figure3.Figure 3
Numerical procedure for nonlinear sloshing simulation.


## 4. Numerical Results and Discussion

A code is developed following the above numerical formulation for computing sloshing response. The liquid sloshing response inside a 2.0 m wide rectangular rigid container with liquid filled to a depth of 1 m is simulated. Aspect ratio of 0.5 is maintained. In the numerical simulations, 40 nodes along thex-direction and 20 nodes along the z-direction are taken.


### 4.1. Free Vibration Analysis

To validate the code for stiffness matrix formulation, a free vibration problem is solved first. A mass matrixM (shown in (24)) for the free surface of the liquid is formed:
(24)M=1g∫ΓsNTNdΓs.
If ωn denotes the nth natural slosh frequency of the coupled system and {Ψn} the corresponding mode shape, the free vibration problem to be solved is
(25)(K-ωn2M){Ψn}=0.The natural slosh frequencies obtained from (25) are compared with Faltinsen’s analytical solution [4]. For a rectangular tank the order of natural sloshing frequency is [4],
(26)ωn=gkntanh(knh)(n=1,2,3,…),
where kn is the wave number, given by kn=nπ/L, L is the length of the container, and h is the water depth. Table 1 shows the slosh frequencies in rad/s obtained for the present tank using finite element method and the above analytical formula. Both the results are in good match. Figure 4 shows the first four mode shapes obtained.Table 1
Slosh frequencies compared with analytical solution.
Mode no.
Present(rad/s)
Equation (26) (rad/s)
Error %
1
3.7607
3.7594
0.0353
2
5.5456
5.5411
0.0801
3
6.8092
6.7986
0.1553
4
7.8715
7.8510
0.2620
5
8.8126
8.7777
0.3983First four mode shapes of sloshing.
(a)
(b)
(c)
(d)


### 4.2.  Sloshing Response under Harmonic Excitation

In this section, sloshing response is simulated, when the tank is excited with harmonic motion. The tank is assumed to be subjected to the following forced harmonic horizontal motion:(27)Xt(t)=ahcos(ωht),
where ah is horizontal forcing amplitude, t is time, and ωh is the angular frequency of the forced horizontal motion. Equation (27) gives excitation velocity as -ahωhsin(ωht), which leads to a zero-free surface velocity potential initial condition from (7). In the present numerical simulation, 40 nodes along the x-direction and 20 nodes along the z-direction are taken, and a time step of 0.003 s is adopted. The sloshing response is evaluated for various excitation frequencies which are closer, away, and equal to the fundamental sloshing frequency. The free surface behaviour is examined for smaller and steeper wave according to Frandsen [10]. The simulation cases considered are shown in Table 2. The results obtained are compared with numerical results of Frandsen [10].Table 2
Simulation cases for the sloshing response.
Case
Frequencyωh (rad/s)
Amplitudeah (m)
Smaller wave
Steeper wave
1
0.7ω1
0.0052
0.0520
2
0.9ω1
0.0039
0.0393
3
ω1
9.9 * 10−4
0.0097
4
1.3ω1
0.0015
0.0301Figures5(a) and 5(b) show the free-surface elevation at the left wall for case 1 for smaller and large horizontal forcing amplitude. The external excitation frequency is in off resonance region; the forcing frequency is smaller than the first fundamental sloshing frequency. The time histories of the sloshing response are nondimensionalised with the first natural sloshing frequency. The results obtained are compared with Frandsen and found to be in good agreement. Figures 5(c) and 5(d) show the associated phase plane plots. These phase plane plots display linear behaviour of the free surface.Time history response and associated phase plane plot for case 1.
(a)
(b)
(c)
(d)Figures6(a) and 6(b) show the slosh response for case 2 for small and large forcing amplitude, respectively. In this case the excitation frequency is closer to the first slosh natural frequency and as expected a beating phenomenon is observed. For small forcing amplitude the response is compared with linear solution and both the results are in good match the response is symmetric in this case displaying a linear behaviour. For large forcing amplitude, the system is in nonlinear region; due to nonlinearities an asymmetric beating phenomenon is observed.Time history response for case 2.
(a)
(b)Figures7(a) and 7(b) show the slosh response at the left wall for case 3 for small and large forcing amplitude, respectively, and the corresponding phase plane plots. In this case, the external excitation frequency is equal to the fundamental slosh frequency and as expected resonance takes place. The results obtained are in excellent agreement with the Frandsen results. For small amplitude case the response is almost symmetric, but for large amplitude case the response is not symmetric because of nonlinear effects. The phase plane plots show clearly the difference between the small forcing amplitude and large forcing amplitude. A moving mesh generated at different type steps for large forcing amplitude case is shown in Figure 8. (Animation of simulation in this case can be seen in http://www.youtube.com/watch?v=LlwUOWMmVtc.) It helps in understanding the sloshing flow patterns.Time history response and associated phase plane plots for case 3.
(a)
(b)
(c)
(d)Generated mesh in case 3 for large forcing amplitude at various time steps.
(a)
(b)
(c)
(d)
(e)
(f)Figures9(a) and 9(b) show the slosh response at the left wall for case 4. This case is also an off resonance case as case 1 but the forcing frequency is higher than the first natural sloshing frequency. The results obtained are in good agreement with Frandsen.Time history response for case 4.
(a)
(b)


### 4.3. Sloshing Response under Random Excitation

In this section sloshing response is simulated for random excitations. For simulating random sloshing, first a random excitation time history is needed. The required random excitation is generated using Bretschneider spectrum:(28)Sη=5Hs216ωp(ωpω)5exp[-54(ωpω)4],
where Hs is the significant wave height and ωp is the peak frequency. Since the higher frequencies have no influence on the sloshing waves, a cut of frequency for random waves is set. In this simulation the cut of frequency is taken as five times the natural slosh frequency. Based on this spectrum, the random waves are generated which are given as the base excitation to the container. The displacement of the random wave can be obtained by linear superposition of a series of harmonic waves with random phase as a time function:
(29)x(t)=∑i=1NωAisin(ωit+ϕi),
where x(t) denotes a random horizontal oscillation that the container is subjected to, t is time, ωi is the frequency of ith linear wave, and Nω is the number of all the linear harmonic waves. The frequency ωi ranges from 0 to π/dt. Ai and ϕi are the wave amplitude and phase of each linear harmonic wave, respectively. The wave amplitude is determined by the following equation:
(30)Ai=2Sη(ω)Δω,
where Δω is the frequency interval. The phase ϕi is a random variable and uniformly distributed in the interval [0,2π].The specified spectrum of oscillation withHs=0.01 h and peak frequency ωp=ω1 is shown in Figure 10. To generate random wave from the shown spectrum, Nω is taken as 512, 1024 data points are taken, and a time step of 0.02 s is adopted. The horizontal oscillation of the container corresponding to the spectrum is shown in Figure 11(a). Figure 11(b) shows the corresponding time history of the slosh response at the left wall of the container.Figure 10
Spectrum of oscillation of the tank withHs=0.01 h, ωp=ω1.(a) Random horizontal excitation, (b) slosh response at the left wall of the container.
(a)
(b)


## 4.1. Free Vibration Analysis

To validate the code for stiffness matrix formulation, a free vibration problem is solved first. A mass matrixM (shown in (24)) for the free surface of the liquid is formed:
(24)M=1g∫ΓsNTNdΓs.
If ωn denotes the nth natural slosh frequency of the coupled system and {Ψn} the corresponding mode shape, the free vibration problem to be solved is
(25)(K-ωn2M){Ψn}=0.The natural slosh frequencies obtained from (25) are compared with Faltinsen’s analytical solution [4]. For a rectangular tank the order of natural sloshing frequency is [4],
(26)ωn=gkntanh(knh)(n=1,2,3,…),
where kn is the wave number, given by kn=nπ/L, L is the length of the container, and h is the water depth. Table 1 shows the slosh frequencies in rad/s obtained for the present tank using finite element method and the above analytical formula. Both the results are in good match. Figure 4 shows the first four mode shapes obtained.Table 1
Slosh frequencies compared with analytical solution.
Mode no.
Present(rad/s)
Equation (26) (rad/s)
Error %
1
3.7607
3.7594
0.0353
2
5.5456
5.5411
0.0801
3
6.8092
6.7986
0.1553
4
7.8715
7.8510
0.2620
5
8.8126
8.7777
0.3983First four mode shapes of sloshing.
(a)
(b)
(c)
(d)


## 4.2.  Sloshing Response under Harmonic Excitation

In this section, sloshing response is simulated, when the tank is excited with harmonic motion. The tank is assumed to be subjected to the following forced harmonic horizontal motion:(27)Xt(t)=ahcos(ωht),
where ah is horizontal forcing amplitude, t is time, and ωh is the angular frequency of the forced horizontal motion. Equation (27) gives excitation velocity as -ahωhsin(ωht), which leads to a zero-free surface velocity potential initial condition from (7). In the present numerical simulation, 40 nodes along the x-direction and 20 nodes along the z-direction are taken, and a time step of 0.003 s is adopted. The sloshing response is evaluated for various excitation frequencies which are closer, away, and equal to the fundamental sloshing frequency. The free surface behaviour is examined for smaller and steeper wave according to Frandsen [10]. The simulation cases considered are shown in Table 2. The results obtained are compared with numerical results of Frandsen [10].Table 2
Simulation cases for the sloshing response.
Case
Frequencyωh (rad/s)
Amplitudeah (m)
Smaller wave
Steeper wave
1
0.7ω1
0.0052
0.0520
2
0.9ω1
0.0039
0.0393
3
ω1
9.9 * 10−4
0.0097
4
1.3ω1
0.0015
0.0301Figures5(a) and 5(b) show the free-surface elevation at the left wall for case 1 for smaller and large horizontal forcing amplitude. The external excitation frequency is in off resonance region; the forcing frequency is smaller than the first fundamental sloshing frequency. The time histories of the sloshing response are nondimensionalised with the first natural sloshing frequency. The results obtained are compared with Frandsen and found to be in good agreement. Figures 5(c) and 5(d) show the associated phase plane plots. These phase plane plots display linear behaviour of the free surface.Time history response and associated phase plane plot for case 1.
(a)
(b)
(c)
(d)Figures6(a) and 6(b) show the slosh response for case 2 for small and large forcing amplitude, respectively. In this case the excitation frequency is closer to the first slosh natural frequency and as expected a beating phenomenon is observed. For small forcing amplitude the response is compared with linear solution and both the results are in good match the response is symmetric in this case displaying a linear behaviour. For large forcing amplitude, the system is in nonlinear region; due to nonlinearities an asymmetric beating phenomenon is observed.Time history response for case 2.
(a)
(b)Figures7(a) and 7(b) show the slosh response at the left wall for case 3 for small and large forcing amplitude, respectively, and the corresponding phase plane plots. In this case, the external excitation frequency is equal to the fundamental slosh frequency and as expected resonance takes place. The results obtained are in excellent agreement with the Frandsen results. For small amplitude case the response is almost symmetric, but for large amplitude case the response is not symmetric because of nonlinear effects. The phase plane plots show clearly the difference between the small forcing amplitude and large forcing amplitude. A moving mesh generated at different type steps for large forcing amplitude case is shown in Figure 8. (Animation of simulation in this case can be seen in http://www.youtube.com/watch?v=LlwUOWMmVtc.) It helps in understanding the sloshing flow patterns.Time history response and associated phase plane plots for case 3.
(a)
(b)
(c)
(d)Generated mesh in case 3 for large forcing amplitude at various time steps.
(a)
(b)
(c)
(d)
(e)
(f)Figures9(a) and 9(b) show the slosh response at the left wall for case 4. This case is also an off resonance case as case 1 but the forcing frequency is higher than the first natural sloshing frequency. The results obtained are in good agreement with Frandsen.Time history response for case 4.
(a)
(b)


## 4.3. Sloshing Response under Random Excitation

In this section sloshing response is simulated for random excitations. For simulating random sloshing, first a random excitation time history is needed. The required random excitation is generated using Bretschneider spectrum:(28)Sη=5Hs216ωp(ωpω)5exp[-54(ωpω)4],
where Hs is the significant wave height and ωp is the peak frequency. Since the higher frequencies have no influence on the sloshing waves, a cut of frequency for random waves is set. In this simulation the cut of frequency is taken as five times the natural slosh frequency. Based on this spectrum, the random waves are generated which are given as the base excitation to the container. The displacement of the random wave can be obtained by linear superposition of a series of harmonic waves with random phase as a time function:
(29)x(t)=∑i=1NωAisin(ωit+ϕi),
where x(t) denotes a random horizontal oscillation that the container is subjected to, t is time, ωi is the frequency of ith linear wave, and Nω is the number of all the linear harmonic waves. The frequency ωi ranges from 0 to π/dt. Ai and ϕi are the wave amplitude and phase of each linear harmonic wave, respectively. The wave amplitude is determined by the following equation:
(30)Ai=2Sη(ω)Δω,
where Δω is the frequency interval. The phase ϕi is a random variable and uniformly distributed in the interval [0,2π].The specified spectrum of oscillation withHs=0.01 h and peak frequency ωp=ω1 is shown in Figure 10. To generate random wave from the shown spectrum, Nω is taken as 512, 1024 data points are taken, and a time step of 0.02 s is adopted. The horizontal oscillation of the container corresponding to the spectrum is shown in Figure 11(a). Figure 11(b) shows the corresponding time history of the slosh response at the left wall of the container.Figure 10
Spectrum of oscillation of the tank withHs=0.01 h, ωp=ω1.(a) Random horizontal excitation, (b) slosh response at the left wall of the container.
(a)
(b)


## 5. Conclusion

Sloshing response of liquid in 2D fixed and forced tanks is investigated numerically considering fully nonlinear equations. A mixed Eulerian-Lagrangian nonlinear finite element numerical model has been developed based on potential flow theory. Free-surface sloshing response is simulated under regular harmonic base excitations for small and steep waves and then the formulation is extended to random excitations. In the simulations, the tank is assumed to be rigid with aspect ratio ofh/L=0.5. For accurate velocity computation, the velocity field is interpolated from the velocity potential according to least square method. The fourth-order Runge-Kutta method is employed to advance the solution in time domain and a regridding technique based on cubic spline is employed to avoid numerical instabilities. The test cases considered here are similar to the test cases by Frandsen for regular harmonic excitations. The results obtained are in perfect match with Frandsen results. In the present simulation, no regridding is required for the simulation of small amplitude waves. In the case of steep waves, for long time simulation, regridding should be carried out for every few time steps.



---

*Source: 571528-2013-02-27.xml*