SPH and ALE formulations for sloshing tank analysis

Design of fuel tanks requires the knowledge of hydrodynamic pressure distribution on the structure. These can be very useful for engineers and designers to define appropriate material properties and shell thickness of the structure to be resistant under sloshing or hydrodynamic loading. Data presented in current tank design codes as Eurocode, are based on simplified assumptions for the geometry and material tank properties. For complex material data and complex tank geometry, numerical simulations need to performed in order to reduce experimental tests that are costly and take longer time to setup. Different formulations have been used for sloshing tank analysis, including ALE (Arbitrary Lagrangian Eulerian) and SPH (Smooth Particle Hydrodynamic). The ALE formulation uses a moving mesh with a mesh velocity defined trough the structure motion. In this paper the mathematical and numerical implementation of the FEM and SPH formulations for sloshing problem are described. From different simulations, it has been observed that for the SPH method to provide similar results as ALE formulation, the SPH meshing, or SPH particle spacing needs to be finer than ALE mesh. To validate the statement, we perform a simulation of a sloshing analysis inside a partially filled tank. For this simple, the particle spacing of SPH method needs to be at least two times finer than ALE mesh. A contact algorithm is performed at the fluid structure interface and SPH particles. In the paper the efficiency and usefulness of two methods, often used in numerical simulations, are compared.


INTRODUCTION
In computational mechanics there is a growing interest in developing meshless methods and particle methods as alternatives to traditional FEM (Finite Element Method).Among the various meshfree and particle methods, Smoothed Particle Hydrodynamics (SPH) is the longest established and is approaching a mature stage.SPH is a Lagrangian meshless method in which the problem to be solved is discretized using particles that are free to move rather than element tied by classical mesh connectivity.SPH method has been extensively used for high impact velocity applications, in aerospace and defense industry for problems where classical FEM methods fail due to high meshes distortion.For small deformation, FEM Lagrangian *Corresponding Author: E-mail: mhamed.souli@univ-lille1.frformulation can solve structure interface and material boundary accurately, the main limitation of the formulation is high mesh distortion for large deformation and moving structure.One of the commonly used approach to solve these problems is the ALE formulation which has been used with success for simulation of fluid structure interaction with large structure motion such as sloshing fuel tank in automotive industry and bird impact in aeronautic industry.It is well known from previous analysis, see Aquelet et al. [1], that the classical FEM Lagrangian method is not suitable for most of the FSI problems due to high mesh distortion in the fluid domain.To overcome difficulties due to large mesh distortion, ALE formulation has been the only alternative to solve fluid structure interaction for engineering problems.For the last decade, SPH and DEM (Discrete Element Method), have been used usefully for engineering problems to simulate high velocity impact problems, high explosive detonation in soil, underwater explosion phenomena, and bird strike in aerospace industry, see Han et al. [2] for details description of DEM method.SPH is a mesh free Lagrangian description of motion that can provide many advantages in fluid mechanics and also for modelling large deformation in solid mechanics.For some applications, including underwater explosion and hydrodynamic impact on deformable structures, engineers have switched from ALE to SPH method to reduce CPU time and save memory allocation.Unlike ALE method, and because of the absence of the mesh, SPH method suffers from a lack of consistency than can lead to poor accuracy, as described in Randles et al. [14] and Vignjevic et al. [16].
In this paper, devoted to ALE and SPH formulations for sloshing tank problems, the mathematical and numerical implementation of FEM and SPH formulations are described.From different simulations, it has been observed that for the SPH method to provide similar results as FEM Lagrangian formulations, the SPH meshing, or SPH spacing particles needs to be finer than the ALE mesh, see Messahel et al. [15] for underwater explosion problem.To validate the statement on fluid structure interaction problems, we perform a simulation of a sloshing problem.In the simulation, the particle spacing of SPH method needs to be at least two times finer than ALE mesh.A contact algorithm is performed at the fluid structure interface for both SPH and ALE formulations.
In Section 2, the governing equations of ALE formulation are described.In this section, we discuss mesh motion as well as advection algorithms used to solve mass, momentum and energy conservation in ALE formulation.Section 3 describes the SPH formulation, unlike FEM formulation which based of the Galerkin approach, SPH is a collocation method.The last section is devoted to numerical simulation of fluid sloshing inside a moving rigid structure tank, using both FEM and SPH methods.To get comparable results between FEM and SPH, the particle spacing of SPH method needs to be at least two times finer than FEM mesh.

ALE FORMULATION
A brief description of the FEM formulation used in this paper is presented, additional details can be provided in Benson [4].To solve fluid structure interaction problems, a Lagrangian formulation is performed for the structure and an ALE formulation for the fluid material.In general ALE description, an arbitrary referential coordinate is introduced in addition to the Lagrangian and Eulerian coordinates.The material derivative with respect to the reference coordinate can be described in equation (2.1).Thus substituting the relationship between material time derivative and the reference configuration time derivative leads to the ALE equations, (2.1) where X i is the Lagrangian coordinate, x i the Eulerian coordinate, w i is the relative velocity.Let denote by v the velocity of the material and by u the velocity of the mesh.In order to simplify the equations we introduce the relative velocity w = v − u.Thus the governing equations for the ALE formulation are given by the following conservation equations: (i) Mass equation. ( (ii) Momentum equation.
σ ij is the stress tensor defined by σ = −p + τ, where τ is the shear stress from the constitutive model, and p the pressure.The volumetric compressive stress p is computed though an equation of state, and the shear stress from material constitutive law.
(iii) Energy equation. ( Note that the Eulerian equations commonly used in fluid mechanics by the CFD community, are derived by assuming that the velocity of the reference configuration is zero, u = 0, and that the relative velocity between the material and the reference configuration is therefore the material velocity, w = v.The term in the relative velocity in (2.3) and (2.4) is usually referred to as the advective term, and accounts for the transport of the material past the mesh.It is the additional term in the equations that makes solving the ALE equations much more difficult numerically than the Lagrangian equations, where the relative velocity is zero.There are two ways to implement the ALE equations, and they correspond to the two approaches taken in implementing the Eulerian viewpoint in fluid mechanics.The first way solves the fully coupled equations for computational fluid mechanics; this approach used by different authors can handle only a single material in an element as described for example in Benson [4].The alternative approach is referred to as an operator split in the literature, where the calculation, for each time step is divided into two phases.First a Lagrangian phase is performed, in which the mesh moves with the material, in this phase the changes in velocity and internal energy due to the internal and external forces are calculated.The equilibrium equations are: (2.5) (2.6)In the Lagrangian phase, mass is automatically conserved, since no material flows across element boundaries.
In the second phase, the advection phase, transport of mass, energy and momentum across element boundaries are computed; this may be thought of as remapping the displaced mesh at the Lagrangian phase back to its original for Eulerian formulation or arbitrary location for ALE formulation using smoothing algorithms.From a discretization point of view of (2.5) and (2.6), one point integration is used for efficiency and to eliminate locking, Belytschko et al. [3], zero energy modes are controlled with an hourglass viscosity.A shock viscosity, with linear and quadratic terms developed by Von-Neumann and Richtmeyer in earlier fifties is used to resolve the shock wave.The resolution is advanced in time with the central difference method, which provides a second order accuracy for time integration.
For each node, the velocity and displacement are updated as follows: (2.7) where F int is the internal vector force and F ext the external vector force associated with body forces, coupling forces, and pressure boundary conditions, M is a diagonal lumped mass matrix.For each element of the mesh, the internal force is computed as follows: B is the gradient matrix and Nelem the number of elements.The time step size, Δt, is limited by the Courant stability condition [4], which may be expressed as: where l is the characteristic length of the element, and c the sound speed of the element material.For a solid material, the speed of sound is defined as: where ρ is the material density, K is the module of compressibility.

MOVING MESH ALGORITHM
The ALE algorithm used in the paper allows the fluid nodes to move in order to maintain the integrity of the mesh.As the fluid impacts the plate, the fluid mesh moves with a mesh velocity that is different from fluid particle velocity.The choice of the mesh velocity constitutes one of the major problems with the ALE description.Different techniques have been developed for updating the fluid mesh domain.For problems defined in simple domains, the mesh velocity can be deduced through a uniform or non uniform distribution of the nodes along straight lines ending at the moving boundaries.This technique has been used for different applications including water wave problems.For general computational domains, the mesh velocity is computed through partial differential equations, with appropriate boundary conditions.For sloshing problem where the tank is moving with a applied velocity that is time dependent, the fluid mesh moves as a rigid mesh following the tank.This new ALE feature allows the mesh to stay regular, and the time step, which can be affected by mesh distortion, to be stable.In other words, there is only mesh motion and no mesh distortion due to the ALE formulation.This method is very useful for moving or rotating tanks, where the fluid mesh will move and rotate with the tank without undergoing any mesh deformation.
The ALE algorithm used in the paper allows the fluid mesh to follow the movement of the structure.The integrity of mesh structure is maintained.As the structure impacts the rigid plate and then moves and rotates, the fluid mesh moves as a rigid mesh in the coordinate system attached to the structure.This ALE algorithm can be applied to several problems in moving structure that are rigid or undergo small deformations.

ADVECTION PHASE
In the second phase, the transport of mass, momentum and internal energy across the element boundaries is computed.This phase may be considered as a 're-mapping' phase.The displaced mesh from the Lagrangian phase is remapped into the initial mesh for an Eulerian formulation.To illustrate the advection phase, we consider in Figure 1.1, a simple problem with 2 different materials, one with high pressure and the second a lower pressure.During the Lagrangian phase, material with high pressure expands, and the mesh moves with the material.Since we are using Eulerian formulation, the mesh is mapped to its initial configuration, in the advection phase, material volume called flux is moving from element to element, but we keep separate materials in the same element, using interface tracking between materials inside an element.Conservation properties are performed during the Lagrangian phase; stress computation, boundary conditions, contact forces are computed.The advection phase can be seen as a remapping phase from a deformable mesh to initial mesh for an Eulerian formulation, or to an arbitrary mesh for general ALE formulation.In the advection phase, volume flux of material through element boundary needs to be computed.
Once the flux on element faces of the mesh is computed, all state variables are updated according to the following algorithm, using a finite volume algorithm (2.10), (2.10) where the superscripts '−' and '+' refer to the solution values before and after the transport.Values that are subscripted by j refer to the boundaries faces of the element, through which the material flows, and the Flux j are the volume fluxes transported through adjacent elements.The flux is positive if the element receives material and negative if the element looses material.
The ALE multi-material method is attractive for solving a broad range of non-linear problems in fluid and solid mechanics, because it allows arbitrary large deformations and enables free surfaces to evolve.The advection phase of the method can be easily implemented in an explicit Lagrangian finite element code.

SPH FORMULATION
The SPH method developed originally for solving astrophysics problem has been extended to solid mechanics by Libersky et al. [8] to model problems involving large deformation including high impact velocity.SPH method provides many advantages in modeling severe deformation as compared to classical FEM formulation which suffers from high mesh distortion.The method was first introduced by Lucy [9] and Monaghan et al. [11], for gas dynamic problems and for problems where the main concern is a set of discrete physical particles than the continuum media.The method was extended to solve high impact velocity in solid mechanics, CFD applications governed by Navier-Stokes equations, and fluid structure interaction problems.It is well known from previous papers, Vila [17] that SPH method suffers from lack of consistency that can lead to poor accuracy of motion approximation.
A detailed overview of the SPH method is developed by Liu et al. [10], where the two steps for representing a continuous function, using integral interpolation and kernel approximation are given by (3.1) and (3.2), (3.1)where the Dirac function satisfies: The approximation of the integral function (3.1) is based on the kernel approximation W, that approximates the Dirac function based on the smoothing length h, Unlike FEM, where weak form Galerkin method with integration over mesh volume, is common practice to obtain discrete set of equations, in SPH since there is no mesh, ∑ ) / ( h d θ collocation based method is used.In collocation method the discrete equations are obtained by enforcing equilibrium equations, mass, linear momentum and energy, at each particle.In SPH method, the following equations are solved for each particle (3.3), (3.3)where ρ i , v i , e i are density, velocity and internal energy of particle, σ i , m i are Cauchy stress and particle mass.
A ij is the gradient of the kernel function defined by (3.4) It had be shown that convergence and stability of the SPH methods depends on the distribution of particles in the domain.
In order to treat problem involving discontinuities in the flow variables such as shock waves an additional dissipative term is added as an additional pressure term.This artificial viscosity should be acting in the shock layer and neglected outside the shock layer.In this simulation a pseudo-artificial pressure term μ ij derived by Monoghan et al. [12] is used.This term is based on the classical Von Neumann artificial viscosity and is readapted to the SPH formulation as follow, Where are respectively the average density and speed of sound, ∈ is small perturbation that is added to momentum and energy equations to avoid singularities, finally α and β in equation (5.1) are respectively the linear and quadratic coefficients.

CONSTITUTIVE MATERIAL MODELS FOR WATER
In this paper a Newtonian fluid constitutive material law is used water, see Where C is the bulk speed of sound, where ρ 0 and ρ are the initial and current densities.The coefficient s is the linear Hugoniot slope coefficient of the shock velocity particle velocity (U s − U p ) curve, equation (4.2), U s is the shock wave velocity, U p is the particle velocity.This equation of state requires fluid specific coefficient S, which is obtained through shock experiment by curve fitting of the U s − U p relationship.Shock experiments on fluids and solids provide a relation between the shock speed U s and the particle velocity behind the shock, U p along the locus of shocked states.
An important phenomenon that arises during hydrodynamic impact is the formation of shock, mathematically equations (3.1)-(3.3)develop a shock, which lead to non continuous solution and the problem is well posed only if the shock conditions are satisfied.These conditions called Rankine Hugoniot conditions describe the relationship between the states on both sides of the shock for conservation of mass, momentum and energy across the shock, and are derived by enforcing the conservation laws in integral form over a control volume that includes the shock.In the absence of numerical viscosity, high non physical oscillations are generated in the immediate vicinity of the shock.

FLUID-STRUCTURE CONTACT ALGORITHM
For SPH and ALE simulations, a penalty type contact is used to model the interaction between the fluid and the plate.In computational mechanics, contact algorithms have been extensively studied by several authors.Details on contact algorithms can be found in Belyshko et al. [3].Classical implicit and explicit coupling are described in detail in Longatte et al. [12], where hydrodynamic forces from the fluid solver are passed to the structure solver for stress and displacement computation.In this paper, a coupling method based on penalty contact algorithm is used.In penalty based contact, a contact force is computed proportional to the penetration depth, the amount the constraint is violated, and a numerical stiffness value.In an explicit FEM method, contact algorithms compute interface forces due to impact of the structure on the fluid, these forces are applied to the fluid and structure nodes in contact in order to prevent a node from passing through contact interface.In contact one surface is designated as a slave surface, and the second as a master surface.The nodes lying on both surfaces are also called slave and master nodes respectively.The penalty method imposes a resisting force to the slave node, proportional to its penetration through the master segment, as shown in Figure 4.1 describing the contact process.This force is applied to both slave and nodes of the master segment in opposite directions to satisfy equilibrium.Penalty coupling behaves like a spring system and penalty forces are calculated proportionally to the penetration depth and spring stiffness.The head of the spring is attached to the structure or slave node and the tail of the spring is attached to the master node within a fluid element that is intercepted by the structure, as illustrated in Figure 4.1.
Similarly to penalty contact algorithm, the coupling force is described by (5.1): where k represents the spring stiffness, and d the penetration.The force F in Figure 4.1 is applied to both master and slave nodes in opposite directions to satisfy force equilibrium at the interface coupling, and thus the coupling is consistent with the fluid-structure interface conditions namely the action-reaction principle.
The main difficulty in the contact algorithms comes from the evaluation of the stiffness coefficient k in Eq. (5.2).The stiffness value is problem dependent, a good value for the stiffness should reduce energy interface in order to satisfy total energy conservation, and prevent fluid leakage through the structure.The value of the stiffness k is still a research topic for explicit contact-impact algorithms in structural mechanics.In this paper, the stiffness value is similar to the one used in Lagrangian explicit contact algorithms.The value of k is given in term of the bulk modulus K of the fluid element in the coupling containing the slave structure node, the volume V of the fluid element that contains the master fluid node, and the average area A of the structure element connected to the structure node.

ALE MESH SENSITIVITY ANALYSIS FOR SPH METHOD
A detailed finite element model was developed to represent the sloshing problem.Before conducting the simulation, mesh sensitivity tests were performed to compute sloshing frequency of the finite element model for which analytical solution is available in the literature.Three different mesh densities were used for mesh sensitivity tests from 20.000 to 60.0000 hexahedra elements for the fluid mesh.Simulation of the three finite element meshes gives same results with good correlation experimental test results provided by Shao et al. [20].The optimal model of 20.000 elements, shown in Figure 5.2, was taken as reference solution for the ALE finite element simulation.In the simulation, the tank is modeled as rigid material shearing common nodes with the fluid mesh at the fluid structure interface.Dimensions of the rigid tank are presented in Figure 5.2, with a thickness of 1 mm.
For the SPH simulations, three different models have been used, the first model has a number of particles similar to the number of nodes in the ALE model, which consists of The tank is submitted to a horizontal velocity described in the paper by Sho et al. [20] where experimental time history for the height of the waver for 10 seconds is provides.The tank velocity in the horizontal direction is given by: v(t) = 0.032⋅cos(2⋅π ⋅t/T) where T = 1.5 sec is the period of the horizontal velocity motion of the tank.The amplitude of the water height at the tank wall are computed for ALE and SPH simulations and compared to experimental data.Good correlation between ALE and experimental results using identical parameters for the water and tank as shown in Table 3.
The first SPH simulation using 20.000 particles similar to ALE model, did not correlate well with either ALE nor experimental results in term of the amplitude of the wave at the tank interface.To improve SPH results and obtain good correlation with ALE model, finer particle spacing needs to be performed for SPH simulation.SPH refinement can be performed by decreasing particle pacing by a factor from two to four, which can be achieved by increasing the number of SPH particles from 20.000 to 75.00 and 120.00 particles, where both SPH discretizations of two SPH models are shown in Figures 5.3 and 5.4.By refining the SPH model we achieved good correlations between SPH and ALE models in term of the height of the water wave, Figure 5.11 shows time history of height of the water wave at the structure.The price that need to be paid for efficiency of SPH method, is that the SPH method may need larger number of particles to achieve an accuracy comparable with that of a mesh based method.It is well known from design engineers and FEM analysts that experimental tests are costly and take long time to perform.To reduce the number of experimental tests, numerical simulations need to validated and then performed on different design product before setting up a prototype.In order to validate the SPH technique described in the paper, the ALE formulation can be used for validation, since ALE solution is accurate for times where the mesh is deformed but not highly distorted, and has been validated for different applications.
The biggest advantage the SPH method has over ALE methods is that it avoids the heavy tasks of re-meshing.For some complex fluid structure interaction simulations where elements need to be eroded due to failure, the ALE remeshing method may fail, since a new element connectivity needs to be regenerated.SPH method allows failure particles by deactivating failed particles for the particle loop processing.This is a major advantage that SPH method has over classical ALE and classical FEM methods.
To further improve the accuracy of the SPH method for the simulation of free surface and impact problems, efficiency and usefulness of the two methods, often used in numerical simulations, are compared.

CONCLUSION
For structure integrity, several efforts have been made in automotive industry, for modeling sloshing tank analysis and their effect on structures.In automotive and aerospace industries, engineers and FEM analysts move their simulations from mesh based method to SPH method to simulate water impact on deformable structure.We also observed in defense industry where SPH method is recently used for undermine explosion problem and their impact on the surrounding structure.
The biggest advantage the SPH method has over mesh based mesh methods is that it avoids the heavy tasks of re-meshing for hydrodynamic problems or structural problems with large deformation.The price to be paid for efficiency is that the SPH method may need finer resolution to achieve accuracy comparable with that of a mesh based.As a result, SPH simulation can be utilised by using finer particle spacing for applications where mesh based method cannot be used because of remeshing problems due to high mesh distortions.Since the ultimate objective is the design of a safer structure, numerical simulations can be included in shape design optimization with shape optimal design techniques (see Barras et al. [19]), and material optimization (see Gabrys et al. [18]).Once simulations are validated by test results, it can be used as design tool for the improvement of the system structure involved.
Figure 3.1, that represents support of the kernel function.SPH interpolation is given by the following equation (3.2):where is the volume of the particle.The sum in equation (3.2) is over particles in the support domain of the kernel as described in Figure 3.1.

Table 1 .
For pressure response, Mie Gruneisen equation of state is used with parameters defined in Table2, which defines material's volumetric strength and pressure to density ratio.Pressure in Mie Gruneisen equation of state is defined by equation(4.1) (4.1)

Table 1 :
1, Material data for water

Table 3 :
ALE and experimental data of peak wave amplitudes