Sloshing in medium size tanks caused by earthquake studied by SPH

A numerical study of sloshing effects in medium-sized liquid storage tanks subjected to earthquake is briefly presented in the paper. The following issues are considered in the study: the phenomena occurring in tanks during excitation, the amount of pressure change during sloshing, and the effects on the tank itself. The numerical model used in the study is based on the Smoothed Particle Hydrodynamics (SPH) method. The presented method can be used for simulating main nonlinear characteristics of the fluid, such as viscosity, cavitation, wave breaking, and turbulence. The reliability of the model and some of its possibilities are illustrated on a practical example.


Introduction
Liquid storage tanks are all around us, and they are crucial for proper functioning of a modern society.During an earthquake, the mass of fluid moves freely and causes drastic changes in pressure that is exerted on the walls and the roof, which can cause extensive tank damage.In many cases, the collapse of the container can have consequences than exceed by far the value of the container and its contents.For example, the collapse of a water tank, a water tower or a similar structure, leads to disruption of the municipal water supply network and, possibly, to the spread of infections.On the other hand, the collapse of tanks in which flammable materials are stored (oil, gasoline, etc.) can lead to environmental pollution, but also to the spread of extensive uncontrolled fire.The effect of motion of the free liquid surface in a partially filled container, caused by any disturbance, is called sloshing [1].To achieve sloshing, the liquid must have a free surface, where fluid dynamics interacting with the container can significantly change the dynamics of the system.The basic problem of sloshing involves estimation of hydrodynamic pressures and distribution of forces, which have a crucial effect on dynamic stability of the tank [2].During the Darfield earthquake (Mw 7.1) in New Zealand (2010), the roof of the drinking water storage cistern collapsed due to uplift forces caused by sloshing [3].During the Chi-Chi earthquake (Mw 7.4) in Taiwan (1999), reinforced concrete roofs of two water purification tanks situated in the plant partially collapsed due to extensive sloshing [4].Without an appropriate water supply, uncontrolled fires caused much more damage than the earthquake itself in the great San Francisco earthquake (1906, USA), and the Kobe earthquake (1995, Japan).Failure of tanks containing highly flammable products can lead to extensive fires, as registered following the Niigata and Alaska earthquakes (1964), and the earthquake in Turkey (1999) when over 17000 people lost life [5].The list of such disasters is very long.
Sloshing is also a very important phenomenon in spacecraft engineering, nuclear engineering, and naval engineering.Sloshing or moving of cargo, water ballast, or other liquids, can cause a ship to capsize.Sloshing can also affect trucks and aircraft, especially during acceleration or deceleration.To avoid damage due to sloshing, tanker operators avoid filling between 10 % and 70 % of the tank height.It's called the barred fill range.Liquid sloshing may have a significant influence on the response of ground liquid storage tanks during an earthquake, especially when the earthquake frequency coincides with the fundamental resonance frequency of the free surface liquid.Tank sloshing has been studied theoretically, experimentally, and numerically over the past several decades.The first pioneering activities in this area were the works of Hough [6] who tried to solve this problem theoretically.Housner developed a mechanical analogue model [7] which is one of the earliest and most used tools for understanding the fundamentals of sloshing in rectangular and cylindrical liquid storage tanks subjected to horizontal ground acceleration.The effects of various parameters on sloshing action in tanks have so far been experimentally investigated: fill level, type of excitation, tank shape, and internal structure of the tank.However, most of these experiments were performed on scaled models rather than on real liquid tanks.Due to limitations of theoretical and experimental studies, a number of numerical models have been developed and used for the simulation of liquid sloshing.For the problems with free-surface flow, moving boundaries and large geometry deformations, the meshfree methods show great advantages compared to mesh-based methods.One such method is the Smoothed Particle Hydrodynamics (SPH), originally developed for astrophysics problems by Lucy [8] and simultaneously by Gingold and Monaghan [9].In recent times, the SPH method has been used in many fields of science, e.g. for multi-phase flow problems [10], large deformation problems of structure [11], underwater explosions [12], and in other fields.Due to its Lagrangian nature, the SPH method appears to be effective in solving liquid sloshing problems.Many researchers studied these problems in twodimensional or three-dimensional space [2].For planar (2D) problems, Delorme et al. [14] analysed a two-dimensional rectangular tank in rolling motion and the characteristics of impact pressure.Marsh et al. [15] analysed liquid sloshing in different tank shapes.Colagrossi et al. [16] and Gotoh et al. [17,18] analysed the violent sloshing phenomenon using the improved SPH method.In 3D, Vorobyev et al. [19] analysed centralized sloshing in a cylindrical tank, with and without obstacles inside.Iglesias et al [20] investigated, numerically and experimentally, the sloshing problem in rectangular tanks exposed to rolling motion.Celebi and Akyildiz [21] studied liquid sloshing in a moving partially filled rectangular tank.Rafiee et al. [22] studied sloshing under a swaying excitation in a rectangular tank, whose width is much smaller than the length.Chen at al. [23] investigated numerically and experimentally the sloshing behaviour of rectangular and cylindrical liquid tanks.Vijay et al. [24] studied effects of sloshing in overhead liquid storage tanks, and Liu and Lin [25] studied three-dimensional (3D) non-linear sloshing with broken free surfaces.It can be seen from literature review that the SPH method can be used to realistically simulate sloshing in tanks in 2D and 3D.Hence, more research is needed to investigate the effects of the internal tank structures and multi-degrees of freedom tank base motion on liquid sloshing.Also, it has been noticed that most numerical and experimental studies of sloshing are focused on harmonic motions, while very few of them deal with real accelerograms.The SPH method is used in this paper to simulate a liquid (water) moving in medium size storage tanks during an Sloshing in medium size tanks caused by earthquake studied by SPH earthquake.The peak acceleration of each earthquake is set to 0.6g, which is the peak ground acceleration for a 10.000 year return period in Croatia.All accelerograms are scaled according to this value for easier comparison of results.The main aim of this study was to investigate the behaviour of a liquid (sloshing displacements, dynamic pressures, forces and kinetic energy distribution) under different types of ground motion with different durations and periods.

Introduction
Fluid flow is the result of an action generated by externally applied forces.Fluid properties (e.g.density and viscosity) affect fluid flow and are a function of other thermodynamic variables (e.g.pressure and temperature).In many flows, the effects of compressibility and viscosity can be neglected.In this paper, the fluid is considered to be inviscid and weakly compressible.

Governing equations
Fundamental governing equations of fluid dynamics are three fundamental physical principles of conservation: conservation of momentum conservation of mass conservation of energy These laws considered together are better known as the Navier-Stokes equations for incompressible fluid, and can be written as follows [26]: where v α represents the vector of velocity (with α being the coordinate direction -x, y, z), ρ is the mass-density, R α is the sum of mass-volume forces acting on the fluid (also for each direction), P is the pressure and μ is the dynamic viscosity of the fluid.The Navier-Stokes equations govern the motion of fluids and when the fluid flow is complex (including turbulence, cavitation, etc.), then the Navier-Stokes equations can be solved only numerically.If only the gravitational mass force acts on the fluid, the equation (1) can be rewritten in the following form: (2) which directly represents the momentum conservation equation.In the above equation, g is the gravitational acceleration and Θ is the viscous force term.If the fluid is considered inviscid, the so called artificial viscosity (proposed by Gingold and Monaghan [9]) is most frequently used as the viscous force term.The real viscosity of the fluid can be modelled using laminar viscosity or full viscosity via the SPS turbulence model [26,27].According to the continuity equation, changes in fluid density can be calculated with: (3) which represents the mass conservation equation in the Eulerian form.Finally, the energy conservation equation, which represents the first law of thermodynamics, can be written in the following form: (4) where e is the thermal energy and σ α denotes the total stress tensor [26].

SPH method
The smoothed particle hydrodynamics (SPH) is a numerical method used for solving Navier-Stokes equations and for fluid flow simulations.Generally, it is a meshfree Lagrangian particle method, in which the considered volume of fluid is divided into a set of particles.Each particle carries physical determinants such as position, density, mass, pressure, etc.Only a brief description of the method is given in this paper.A more detailed account of the method can be found in papers [26][27][28][29][30].
The fundamental principle of the SPH method is the so called integral interpolation of any quantity function A(r) using the kernel function W(r'-r,h) that approximates a delta function via a compact domain Ω (Figure 1): (5) where r is the radius vector of a considered point and r' is the radius vector of any point in Ω within the smoothing length h that controls the smoothness or roughness of the kernel.This kernel approximation (5) can be discretized using particle approximation of the function A at the particle (interpolation point) "i" with values of the same function in adjacent points "j": (6) Marina Sunara Kusić, Jure Radnić, Nikola Grgić, Alen Harapin The SPH method approximates the particle density using the so called "continuity density" approximation, i.e. using the continuity equation ( 3) in discrete form: (7) via the sum of densities of all neighbouring particles within the region of the kernel function.The total number of neighbouring particles is denoted as N, while the mass and the density of particle "j" are denoted as m j and ρ j , respectively.Also, direct application of the SPH particle approximation to the momentum equation yields the following equation: (8) where Θ ij represents the viscosity term between particles "i" and "j", and can be found in [26,27,29].The selection of the smoothing kernel function is of particular significance, especially in terms of computational accuracy and stability of the method [2,26].One of the most popular kernel functions, which is used in this study, is the cubic function [26][27][28] defined by Monaghan: 1 ≤ q < 2 q = (r'-r)/2h (9) where q is the relative distance between two particles.The SPH method assumes that particle masses and massdensities are known for all particles before the start of the method.Particle mass is a user defined constant, but mass-density is a continuous field of the fluid that must be computed [26][27][28].
It is therefore necessary to include another equation, the socalled equation of state, to link pressure and density.The Tait's equation of state (for water) is adopted in this paper for computing the pressure [1]: (10) where ρ 0 is the reference density (equalling to ρ 0 =1000 kgm -3 for water) and B is the parameter that determines the speed of sound.The thermal energy of each particle is calculated according to Monaghan [31].The standard Predictor-Corrector scheme is adopted [26][27][28] as solver algorithm.Also, the Shepard filter is used for density reinitialisation after each 30 time steps [26,27].The original SPH software was taken from [32] and upgraded by adding the possibility of assigning an arbitrary type of external excitation.The original model is based on the now legendary papers of Lucy [8] and Gingold and Monaghan [9], and has been developed to become a mighty engineering tool that can simulate various phenomena in fluids, such as: viscosity, turbulence, cavitation, sloshing, surface tension, wave breaking, sliding objects, wave impact on structure, etc. [26,27].Also, the SPH method is particularly suited for problems involving large fluid deformations, and highly nonlinear or violent free-surface flows [27].The influence of boundary conditions, multiphase modelling, floating bodies, and interaction between the fluid and moveable/ flexible structure, are still the areas of research within the SPH method.

Boundary treatment
In hydrodynamic problems, two main types of boundaries can be distinguished: free surface boundary and solid wall boundary.Because of the Lagrangian nature of the SPH, there is no need for special treatment of the free surface, and kinematic conditions of the free surface can be considered as automatically satisfied.Dealing with the solid wall boundary is still an open issue in the SPH method [28].Boundary conditions do not appear naturally in the SPH formalism.When a particle approaches a solid boundary, its support domain is truncated, which causes calculation errors and affects accuracy of numerical simulations.Solid wall boundary treatments in SPH method are based on creation of several virtual particles, which simulate the boundary of the computer domain.Three different types of approach can generally be distinguished: the repulsive particles approach (proposed by Gingold & Monaghan [9]), the ghost particles approach (Libersky and Petscheck [33]) and the dynamic particles approach (Crespo et al. [28]).The dynamic particles approach is used in this study.In this approach, the Sloshing in medium size tanks caused by earthquake studied by SPH boundary particles satisfy the same equations as the fluid particles (momentum equation, continuity equation, energy equation, equations of state), but they do not move, i.e. they remain fixed in position (fixed boundaries) or move according to some externally imposed function (moving objects like wavemakers) The most significant advantage of this approach is computational simplicity, since no special considerations are necessary for boundary particles [26][27][28].The friction between fluid (liquid) and walls is not considered because boundary particles only exert a normal force on the fluid.In summary, the issue of boundary conditions is still a point of interest of many researchers.

Numerical tests 3.1. General
The problem of sloshing in medium sized tanks exposed to earthquake forces is investigated in this paper.Such tanks are a reliable choice for long term storage of almost any liquid, ranging from chemicals and fuel to drinking water.These tanks must also exhibit high resistance, which is essential for longterm storage.Dimensions of these tanks vary from small ones for family houses (up to 2.000 l), medium ones for medium size and large buildings (from 2.000 l to 20.000 l), and large ones (over 20.000 l), for all other purposes.In large cities, medium size tanks can often be found on top of the buildings.The behaviour of rectangular and cylindrical tanks is considered in this numerical test (Figure 2).The walls of the tanks are referred to as side walls and the roof as the top wall.The height of the tanks is 3.0 m, and the base surfaces are calculated under the condition that both tanks have the same volume.These tanks are of ordinary dimensions for medium size tanks [34].The tanks are filled with water up to 2 m and 2.85 m, which corresponds to 67 % and 95 % of the tank height.These levels have been chosen because the level of 67 % offers great fluid displacement, but small force on the roof, while level 95 % offers small displacement but great force on the roof [35].The density of ρ 0 =1000 kg/m 3 has been adopted.Using equation of state, the real sound speed c 0 =1430 m/s is corrected with coefficient of the speed of sound, which is 35 (recommended value for water) [31].Also, artificial viscosity is used, with the coefficient of artificial viscosity amounting to 0.1.Similar test for 2D tanks is given in [35].Marina Sunara Kusić, Jure Radnić, Nikola Grgić, Alen Harapin Sloshing in medium size tanks caused by earthquake studied by SPH analytical results can be the consequence of different viscosity treatments.

Table 1. Correlation coefficients between numerical and analytical results for sloshing displacement
Apart from visual comparisons shown in figures 7 and 8, correlation coefficients for sloshing displacement are given in Table 1 for numerical and analytical results.It can be seen that better agreement between numerical and analytical results has been obtained for rectangular tank, compared to cylindrical tank.The poorest agreement has been obtained for sloshing displacement during the Northridge earthquake, which causes the largest sloshing displacement and significantly nonlinear behaviour of the liquid.The main reason for discrepancies between numerical and analytical results lies in different assumptions (analytical solution has been obtained under the assumption of incompressible, irrotational and inviscid liquid).
The large sloshing displacement can occur when the ground motion frequency is close to natural frequencies of the tank fluid.Natural frequencies of rectangular prismatic tanks depend on the tank Each fluid particle is subjected to constant gravity acceleration in z direction (g = -9.81m/s 2 ) and to horizontal ground acceleration in x direction.Four earthquakes (accelerograms) have been selected for analysis.Two of them are the well-known: Kobe (Japan) and Northridge (USA).The other two are: Ston (Croatia) and Banja Luka (Bosnia and Herzegovina).These accelerograms differ greatly according to duration, mode of action and periods, which was the intention during their selection.All accelerograms are scaled to 0.6g to enable easier comparison of numerical test results.The scaled accelerograms are shown in Figure 3. Figure 4 shows the spectrum of all considered earthquakes, with 0.5 % damping.By analysing the acceleration and displacement spectra, a conclusion can be made in advance about the intensity of sloshing displacements in liquid storage tanks.Usually, the free surface movement and response displacements follow this pattern of behaviour, as clearly shown in results given in Section 3.2.
The initial spacing of particles is 0.05 m and they are organized in a Cartesian grid.In the rectangular tank with the fill level of 67 % of the volume, water was simulated with 139.240 particles, while in the rectangular tank with the fill level of 95 % of the volume, water was simulated with 198417 particles (Figure 5).In the cylindrical tank with the fill level of 67 % of the volume, water was simulated with 151280 particles, and in the cylindrical tank with the fill level of 95 % with 215574 particles (Figure 6).The initial time step was set to 5⋅10 -5 s in order to avoid numerical instabilities.The total computation time was approximately 16 days for each example.The computer used for these calculations is a desktop PC (Intel Core i3-3220 3.3 GHz processor with 8 GB of memory).

Results
The easiest way to observe sloshing intensity is by monitoring temporal variation of displacements on the side walls of containers.Sloshing displacement for the rectangular tank with the fill level of 67 % is shown in Figure 7. Sloshing displacement for the cylindrical tank with the same fill level is shown in Figure 8. Results obtained with the previously described SPH model are compared to the Veletsos analytical model for rectangular and cylindrical tanks [36].A good agreement between numerical and analytical results can be seen in terms of sloshing periods and sloshing intensity, especially in the beginning of the analysis.Later on, the damping in the numerical model calms the liquid, and so the oscillation amplitudes become smaller, which is not the case in the analytical model.This is especially evident in the Northridge and Banja Luka earthquakes.Also, the phase shift visible in comparisons between numerical and Marina Sunara Kusić, Jure Radnić, Nikola Grgić, Alen Harapin geometry and fill depth and can be calculated from linear theories [37]: where g is the gravitational acceleration, L the tank length, D the water depth and n the mode number.A similar formula can be found for cylindrical tanks.As seen is Eq.11, an infinite number of natural frequencies can be calculated, but only the first (fundamental) frequency is significant for engineering applications [37,38].
The fundamental natural period of water in the 67 % filled tank is approximately 2 s [1].According to spectral accelerations and displacements of the earthquakes under study (Figure 4), the most unfavourable effect on sloshing displacement is expected for the Northridge earthquake, which is followed by Kobe, Ston and Banja Luka earthquakes (according to the largest displacement/acceleration at 2 s).The expectations have been confirmed with numerical results (figures 7 and 8).Although all considered accelerograms have the same peak ground acceleration, their response differs due to different predominant periods.That confirms the fact that the peak ground acceleration, which is often the most popular parameter for describing earthquakes, is generally a poor parameter for characterizing the damage potential or peak response in earthquake engineering.The rise of water at tank sides is almost 1.0 m for the Northridge earthquake (Figure 7).An important concern that needs to be emphasized at this point is the rising of sloshing displacements when the predominant period of the earthquake, and the fundamental natural period of the liquid inside tank, are very close to each other.Fundamental natural periods for some different water levels are given in Table 2.
The obtained results are also very analogous to the results given in [35].In that paper the authors analysed sloshing in tanks under seismic excitation, but they used unscaled recorded accelerograms, rather than scaled ones.It can clearly be seen in Figure 9 that numerical results from 2D model [35] and numerical results for 3D model (this paper), obtained for original Northridge (nonscaled) accelerogram, are in better mutual agreement then with the analytical model [36].Discrepancies in 3D and 2D numerical model results are primarily visible in terms of sloshing periods.This may be due to different assumptions regarding the computational domain.In 2D model, the third dimension is assumed to be infinitely extended and in that direction the mass flow is zero.Also, the application of different boundary conditions in SPH can cause different results.Sloshing in medium size tanks caused by earthquake studied by SPH Validation of the original numerical model by similar physical processes (creation of waves by landslides, dam-break problem, wave propagation, and wave-structure interaction) can be found in [39].
The visualization of sloshing in rectangular tank (with dimensions shown in Figure 2), caused by the originally recorded Northridge earthquake, is shown in Figure 10.It can be seen that sloshing displacements along the tank width are uniform and that tank sides do not have a major influence in case of relatively small free-surface oscillations.However, when violent sloshing occurs, fluid behaviour becomes nonlinear (e.g.wave break at t = 3.0 s) and water level changes can be observed in transverse direction.These nonlinearities cannot be simulated with the 2D model and analytical model, and they are the cause of differences between analytical [36] and numerical curves.
The time history of kinetic energy alteration has also been studied in the paper to understand the transfer of kinetic energy from liquid to the tank and vice versa.The primary form of energy transfer is the shearing of the fluid.Energy transfer is greater as the velocity gradients in the flow increase.During sloshing, an increase in velocity gradient can occur due to the wave-to-wall and wave-to-wave interactions.
Many investigations [40,41] show that the wave-to-wave interaction is more effective.Also, it can be observed that the kinetic energy in cylindrical tank is slightly greater compared to the energy in rectangular tank, as can be expected according to sloshing levels.It can therefore be concluded that the tank shape and geometry can direct influence the transfer of kinetic energy.Finally, all kinetic energy induced transforms to thermal energy and vanishes from the system.The impact force on the top wall for the rectangular tank and cylindrical tank with the 95 % fill level is shown in Figure 13 and Figure 14, respectively.The impact force caused by sloshing action exerted on the top wall depends on spectral accelerations and displacements of the earthquakes under study, and on the fundamental natural period of the liquid contained in the tank.For this fill level, the fundamental natural period is about 1.98 s and is somewhat lower compared to the fill level of 67 %.The maximum impact   Marina Sunara Kusić, Jure Radnić, Nikola Grgić, Alen Harapin force on the top wall of the tank was registered for the Northridge earthquake.Kobe and Ston have an almost equal impact force, while Banja Luka has a very small impact on the top wall.This is also expected considering the response spectrum (Figure 4).The maximum impact force on the top wall registered for the Northridge earthquake is close to the weight of the liquid contained in the tank, which is about 250 kN.It is also interesting to note that the impact force on the top wall is greater for the cylindrical tanks compared to rectangular tanks, which is also due to tank geometry (Table 3).The pressure distribution along tank walls at some time points is shown in central plane in Figures 15 and 16 for the Northridge earthquake.Two types of dynamic pressure occur during liquid sloshing: non-impulsive pressures and impulsive pressures [1].
Non-impulsive pressures are dynamic pressures that normally occur in an oscillating liquid.Total pressures (hydrostatic pressures + non-impulsive hydrodynamic pressures) are very similar to standard hydrostatic pressure, as can be seen, for example, at pressure diagrams in Figure 15 and Figure 16 for t=2.3 s and t=3.0 s.Impulsive pressures are pressure pulses caused by impact exerted by a liquid on a solid surface.Impulsive pressures are extremely and highly localized pressures that are commonly associated with traveling waves and hydraulic jumps.These pressures can also be seen in Figure 15 and Figure 16, especially at t=2.7 s.It can be concluded that largest impact pressures occur near the free surface, or at tank-wall intersections.This paper clearly shows that reliable results can be obtained with the SPH method for liquid sloshing in rectangular and cylindrical tanks with rigid walls.However, studies are still needed to investigate the effects of the tank structure deformability, and the multi-degrees-of-freedom tank base motion, on liquid sloshing.

Conclusion
The SPH method is a powerful tool for analysing problems of sloshing in liquid containers.The analysis shows that the method is also quite appropriate for various types of fluid dynamic problems, especially those with a free surface boundary.The tests also show that sloshing displacements and wave amplitudes obtained with the model based on the SPH method are in a good agreement with another numerical model [35].The analysis results also show the dependence of predominant period of the earthquake on sloshing displacement.Sloshing displacements are the greatest when the predominant period of the earthquake, and the fundamental natural period of the liquid inside tank, are very close to each other.Sloshing in medium size tanks caused by earthquake studied by SPH Total pressures and pressure forces on the top wall are shown for various earthquakes.Based on the results of the analysis, it can be concluded that tank roofs can be exposed to very large additional sloshing forces.These forces can, in some cases, exceed the weight of the roof, which must be taken into account during tank design.The presented SPH method formulation is very time consuming, which is the greatest disadvantage of the model.However, this disadvantage is expected to become less pronounced with future developments in computer technology.The realistic physical and numerical description of reality still remains as the primary problem.
It is important to note that the presented analysis covers only the tanks with rigid walls and roofs.Flexible walls and roofs of tanks, depending on the level of their flexibility, can significantly change the image of hydrodynamic pressures on the surfaces.This issue is very interesting and will be explored by the authors in the scope of their oncoming research.

Figure 7 .
Figure 7. Free surface displacement (m) for rectangular tank with fill level of 67 %

Figure 8 .Figure 9 .
Figure 8. Free surface displacement (m) for cylindrical tank with fill level of 67 %

Figure 13 .
Figure 13.Impact force at top wall for rectangular tank

Figure 14 .Figure 15 .
Figure 14.Impact force at top wall for cylindrical tank

Figure 16 .
Figure 16.Pressure distribution along cylindrical-tank walls at some time points during Northridge earthquake