Nonlinear Finite Element Analysis of Sloshing

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.


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 2 Advances in Numerical Analysis 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 twodimensional 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 fournoded 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 (ℎ/ ) of 0.5; ℎ is depth of fluid, 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.

Governing Equations
Consider a rectangular tank fixed in Cartesian coordinate system , which is moving with respect to inertial Cartesian coordinate system 0 0 0 . The origins of this system are at the left end of the tank wall at the free surface and pointing upwards in direction. These two Cartesian systems coincide when the tank is at rest. Figure 1 shows the tank in the moving Cartesian coordinate system along with the prescribed boundary conditions. The tank is assumed to get displaced along -axes, and the tank position at time is given by 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 : 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 (Γ ) we have =0, On the free surface (Γ ) dynamic and kinematic conditions hold; they are given as Rewrite the above equations in the Lagrangian form [16]: where is the free surface elevation measured vertically above still water level, is the horizontal acceleration of the tank, and is the acceleration due to gravity.
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 = 0 and = 0 can be written as: 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).

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 Figure 2. By introducing the finite element shape functions the liquid velocity potential can be approximated as where is the shape function, is the number of nodes in the element, and is nodal velocity potential.
On applying Galerkin residual method to the Laplace equation, we get with matrix defined by Matrix 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.

Velocity Recovery.
To track the free surface (4) need velocities, which are computed from the calculated potential field using 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: wherêis any velocity component (̂or̂), , are the Gauss locations, and 1 , 2 , 3 , 4 are unknowns which need to be evaluated. To evaluate these unknowns, a least square fit is considered between̂and : 4 Advances in Numerical Analysis where is 2 × 2 order Gauss integration points. Then, the four unknown coefficients are determined from four simultaneous equations obtained from Substituting the obtained '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 where is velocity component ( or ). Using the patch recovery technique velocity components and are calculated.

Numerical Time Integration and Free Surface Tracking.
After calculating the velocity at a time step , 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 + ∇ . 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 are known and can be represented in a single variable as where = ( 1 , 2 , . . . , +1 ) , The free surface position and associated velocity potential at the next time step + 1 can be expresses as

Regridding Algorithm.
At the beginning of the numerical simulation, the free surface nodes are uniformly distributed along the -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 is obtained. Then the free surface is divided into segments with the identical arc length. The coordinates of node is denoted as ( , ) ( = 1, 2, . . . , +1) and let the arc length between two successive points and + 1 be . Being a uniform regridding, can be expressed as The coordinates of the nodes ( , ) is a function of the arc length : The cubic spline interpolation is used to calculate the coordinates ( , ) and the velocity potential on the new uniform free surface is also obtained in a similar fashion.

Complete Algorithm for Nonlinear
Sloshing. Including all the steps above, the algorithm for numerical simulation of nonlinear sloshing is as shown in Figure 3.

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 the -direction and 20 nodes along the -direction are taken.

Free Vibration Analysis.
To validate the code for stiffness matrix formulation, a free vibration problem is solved first. A mass matrix (shown in (24)) for the free surface of the liquid is formed: If denotes the th natural slosh frequency of the coupled system and {Ψ } the corresponding mode shape, the free vibration problem to be solved is 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], where is the wave number, given by = / , is the length of the container, and ℎ 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.

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: where ℎ is horizontal forcing amplitude, is time, and ℎ is the angular frequency of the forced horizontal motion. Equation (27) gives excitation velocity as − ℎ ℎ sin( ℎ ), which leads to a zero-free surface velocity potential initial condition from (7). In the present numerical simulation, 40 nodes along the -direction and 20 nodes along thedirection 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].   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. Figures 7(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.
Figures 9(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.

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: where is the significant wave height and 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: where ( ) denotes a random horizontal oscillation that the container is subjected to, is time, is the frequency of th linear wave, and is the number of all the linear harmonic waves. The frequency ranges from 0 to / . and are the wave amplitude and phase of each linear harmonic wave, respectively. The wave amplitude is determined by the following equation: where Δ is the frequency interval. The phase is a random variable and uniformly distributed in the interval [0, 2 ]. The specified spectrum of oscillation with = 0.01 h and peak frequency = 1 is shown in Figure 10. To generate random wave from the shown spectrum, 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.

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 of ℎ/ = 0.5. For accurate velocity computation, the velocity field is interpolated from the velocity potential