ALE and Fluid Structure Interaction for Sloshing Analysis

Liquid containment tanks are, generally, subjected to large deformations under severe earthquake conditions due to coupling forces between tank and the contained liquid. The accurate description of these forces is vital in order to diminish or eliminate the potential risk of tank failure during an earthquake. Yet, analytical formulations derived for the seismic analysis of liquid storage tanks are not capable to capture the complex fluid-structure effects since they include many assumptions and simplifications not only for the behavior of fluid and structure but also for the external excitation. On the other hand, an appropriate numerical method allows us to cope with large displacements of free surface of the fluid, high deformations of the structure and correctly predicts the hydrodynamic forces due to the high-speed impacts of sloshing liquid on a tank wall and roof. For this purpose, a new coupling algorithm based on the penalty formulation of finite element method which computes the coupling forces at the fluidstructure interface is developed in this paper. This algorithm is constructed on a two superimposed mesh systems which are a fixed or moving ALE mesh for fluid and a deformable Lagrangian mesh for structure. The fluid is represented by Navier-Stokes equations and coupled system is solved using an explicit time integration scheme. In order to verify the analysis capability of coupling algorithm for tank problems, numerical method is applied for the analyses of a rigid rectangular tank under harmonic excitation and a flexible cylindrical tank subjected to earthquake motion and numerical results are compared with existing analytical and experimental results. Strong correlation between reference solution and numerical results is obtained in terms of sloshing wave height.


INTRODUCTION
The damage of liquid storage tanks has been very widespread during many earthquakes since the hydrodynamic forces exerted by the contained liquid have not been properly determined and considered in the analysis of such structures. The level of these forces is affected by many features of tank such as tank wall flexibility, support conditions, soil properties and frequency content and amplitude of the earthquake motion. Although earlier seismic analysis the nonlinear dynamic response of unanchored cylindrical liquid storage tanks subjected to strong base excitation considering both large amplitude nonlinear liquid sloshing and fluid-structure interaction. Barton and Parker [16] examined the seismic response of anchored and unanchored cylindrical storage tanks subjected to only one direction of horizontal excitation using finite element method. Both added mass method and 3D finite elements were employed to model fluid inside the tank, however neither material nor geometric nonlinearities were taken into account in the analysis. The application of Arbitrary Lagrangian Eulerian (ALE) algorithm of finite element method for the fluid-structure interaction problems was detailed by Souli et al. [17] and Souli and Zolesio [18]. One very useful reference which describes the ALE technique application for large amplitude sloshing problem in a tank is carried out by Aquelet et al. [19]. The sloshing behaviors of fluid in 3D rigid cylindrical and rectangular tanks subjected to horizontal oscillations were addressed with a numerical and experimental study by Chen et al. [20]. Liu and Lin [21] adopted finite difference method which solves Navier-Stokes equations to study 2-D and 3-D viscous and inviscid liquid sloshing in rectangular tanks and verified the results with the linear analytical solution and experimental data. Mitra et al. [22] used finite element method to solve wave equation to quantify liquid sloshing in partially filled 2D rigid annular, horizontal cylindrical and trapezoidal containers. Razzaghi and Eshghi [23] investigated behavior of anchored and unanchored steel cylindrical tanks under near-field and far-field earthquakes using finite element method.
Several laboratory measurements have been conducted on anchored and unanchored tanks to verify analytical and numerical techniques developed for the seismic analysis of such structures. Kana [24] measured experimentally wall stresses of a cylindrical flexible tank induced by sloshing and inertial loads. Clough [25] studied experimentally on anchored and unanchored board tanks in order to evaluate seismic design procedures. As the continuation of this study, Manos and Clough [26] carried out shaking table tests as well as static tilt tests on scaled tank models. They focused on the higher order out of round response of unanchored tank which induced in addition to rocking cantilever type response. Manos [27] carried out experiments to determine impulsive mode frequencies and base-overturning moments of broad and tall tanks. Tanaka et al. [28] conducted dynamic tests on small and large scale models under earthquake loading in order to investigate elephant foot buckling and side slipping behaviour of cylindrical tanks. This paper presents the description of a fluid-structure coupling algorithm which is developed for the computation of interaction forces at the interface of the two materials. In this algorithm, two superimposed mesh systems are considered, a fixed Eulerian or ALE mesh for the fluid and a deformable Lagrangian mesh for the structure. Fluid flow is represented with Navier-Stokes Equations which are presented in the form of ALE finite element formulation. The free surface of the fluid is tracked using the volume-of-fluid (VOF) method. The momentum advection algorithm is applied to transfer the materials across element boundaries. The penetration of structure node through the fluid at fluid-structure interface is prevented by penalty coupling forces which are determined proportionally to the penetration depth. A second order accurate time marching procedure based on central difference method is used to advance the solution in time. The present numerical method is applied for the analysis of the sloshing in tank problems subjected to external excitations. The accuracy of the numerical method is validated with the existing analytical formulation derived from potential flow theory as well as the experimental data in terms of wave height of the liquid free surface.

THE ANALYSIS OF FLUID-STRUCTURE SYSTEMS BY MEANS OF FINITE ELEMENT METHOD
The main concern in fluid-structure interaction problems is the computation of the fluid forces that act on a rigid or deformable structure. Accurate computation of these forces is essential to provide structural safety. In fact, the application of fluid-structure interaction technology allows to correctly predict the hydrodynamic forces act on a structure due to motion of fluid by solving the hydrodynamic equations and by using an appropriate algorithm to communicate forces between fluid and structure for dynamic equilibrium. The numerical solution algorithms by means of finite element method have been developed using either Lagrangian, Eulerian or Arbitrary Lagrangian Eulerian mesh descriptions to model the fluid motion interacts with structure. In general, the choice of which representation to use depends on the characteristics of the specific problem. Especially, for large deformation problems, the selection of appropriate mesh description is very important to properly simulate the physical phenomenon.
In Lagrangian system, mainly used for structural mechanics, all the nodes of the computational mesh follow the movement of the associated material nodes, since the mesh nodes are embedded in the material nodes. This type of modelling allows materials with history-dependent constitutive relations to be treated, the facile tracking of free surface and the interfaces between the different types of material. However, it has a major draw back, since large distortions of the computational domain are not possible without remeshing frequently. In the Lagrangian formulation, the boundary conditions are easily imposed because the edges of the mesh represent the limits of the physical domain during calculation.
In Eulerian systems, which are most commonly used in fluid dynamic modelling, the computational mesh and associated nodes are fixed in position, but the material nodes are free to move with respect to the Eulerian grid. This enables the ability to handle large distortions in the modelling, but at the cost of flow detail resolution and the accurate definition of the interface.
ALE is a technique that has been developed to combine the positive aspects of the two above procedures, where the computational nodes of the mesh can move with the material nodes as in Lagrangian systems, remain fixed like in Eulerian systems, or move arbitrarily of the material nodes [29]. This enables greater distortions to be handled by the model than if it were a Lagrangian model, with greater accuracy than if it were a Eulerian system.
Tank problems subjected to earthquake motions involve large deformations of free surface as well as moving wall boundaries of the fluid domain, which highly deforms at the same time. For this kind of problems an intermediate formulation which provides mesh of the fluid domain following the boundary motion and preserving the volume shape simultaneously is necessary. The ALE algorithm of finite element method constitutes a very efficient intermediate formulation for the solution of the seismic tank analysis since it copes with large deformations of free surface of the fluid and the structure and correctly predicts the hydrodynamic forces due to coupling of flexible tank and the fluid and the high-speed impacts of sloshing liquid on tank wall and roof.

ALE AND EULERIAN MULTI-MATERIAL FORMULATIONS
In the Eulerian Multi-Material formulation, two or more different materials can be mixed within the same fixed mesh where these materials flow through. A unique material number is assigned to each material and their material boundaries are defined. Material numbers therefore serve a function beyond distinguishing material with different properties. ALE Multi-Material formulation uses the same procedure, but it allows the finite element mesh, which is fixed in the Eulerian Multi-Material formulation, to move independently from the material flow [30]. After moving the mesh a remapping of variables to the new mesh locations is performed with advection methods which are essential for accurate remapping. At each time step, state variables are computed and stored for each material.

ALE REPRESENTATION OF NAVIER-STOKES EQUATIONS
The Arbitrary Lagrangian Eulerian (ALE) approach is based on the arbitrary movement of a reference domain which, additionally to the common material (Lagrangian) domain and spatial (Eulerian) domain, is introduced as a third domain, as detailed in [31]. In this reference domain, which will later on correspond to the finite element mesh, the problem is formulated. The arbitrary movement of the reference frame, accompanied of course by a good "mesh moving algorithm", enables us to rather conveniently deal with moving boundaries, free surfaces, large deformations, and interface contact problems.
The total time derivative of a variable f with respect to a reference coordinate can be described as following: (1) where is the Lagrangian coordinate, is the ALE coordinate, is the particle velocity and is the velocity of the reference coordinate, which will represent the grid velocity for the numerical simulation, and the system of reference will be later the ALE grid. Thus substituting the relationship between the total time derivative and the reference configuration time derivative derives the ALE equations.
Let represent the domain occupied by the fluid particles, and let denote its boundary (Figure1). The equations of mass, momentum and energy conservation for a Newtonian fluid in ALE formulation in the reference domain, are given by: Figure 1 Fluid domain and boundary condition description. (4) where ρ is the density and σ is the total Cauchy stress given by: (5) where p is the pressure and µ is the dynamic viscosity. For compressible flow, Equations (2)-(4) are completed by an equation of state that relies pressure to density and internal energy. The internal energy, e, in Equation (4) is the energy per unit volume.
Equations (2)-(4) are completed with appropriate boundary conditions. The part of the boundary at which the velocity is assumed to be specified is denoted by . The inflow boundary condition is: The traction boundary condition associated with Equation (3) is the conditions on stress components. These conditions are assumed to be imposed on the remaining part of the boundary.
where, is the outward unit normal vector on the boundary. One of the major difficulties in time integration of the ALE Navier-Stokes Equations (2)-(4) is due to the nonlinear term related to the relative velocity ( ) which is usually referred to as the advective term. This term accounts for the transport of the material past the mesh and makes solving the ALE equations much more difficult numerically than the Lagrangian equations, where the relative velocity is zero. For some ALE formulations, the mesh velocity can be solved using a remeshing and smoothing process. In the Eulerian formulation, since, the problem is formulated in the spatial coordinate, the reference frame is fixed: . This assumption eliminates the remeshing and smoothing process, but does not simplify the Equations (2)-(4).
In order to solve ALE formulation of Navier-Stokes Equations (2)-(4), there are two ways, 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. The alternative approach is referred to as an operator split method in the literature ( [32] and [17]) where the calculation, for each time step is divided into two phases. First a Lagrangian phase is performed, using an explicit finite element method, in which the mesh moves with the fluid particle. In the second phase, the displaced mesh from the Lagrangian phase is remapped into the initial mesh for an Eulerian formulation, or an arbitrary distorted mesh for an ALE formulation. In the CFD community, the Lagrangian phase is referred to as a linear Stokes problem. In this phase, the changes in ( ) . r r u r uuuu r r r r v velocity, pressure and internal energy due to external and internal forces are computed. The equilibrium equations for the Lagrangian phase are: Although the continuity equation can be used to obtain the density in a Lagrangian formulation, it is simpler and more accurate to use the integrated form Equation (10) in order to compute the current density ρ [29]: (10) where ρ 0 is the initial density and J is the volumetric strain given by the Jacobian: (11) From a discretization point of view Equation (8) and Equation (9) are computed by using one point integration for efficiency and to eliminate locking [30]. The zero energy modes are controlled with an hourglass viscosity [33]. A shock viscosity, with linear and quadratic terms, is used to resolve the shock wave [34]; a pressure term is added to the pressure in the energy equation (Equation (9)). The resolution is advanced in time with the central difference method, which provides a second order accuracy in time using an explicit method in time.
Central difference method, which is derived from Taylor series expansion, is employed to advance the position of the mesh in time using an explicit method in which the solution progresses without any iteration between consecutive time steps. In order to have a second order accurate scheme in time, the velocity must be staggered with respect to the displacement: The internal nodal forces which are a function of stresses and external forces associated with body forces and boundary conditions (boundary forces, non-reflection boundary forces and contact forces) are used to update the velocity in time: (13) where F int is the internal force vector and F ext is the external force vector and M is the mass matrix diagonalized. For each element of the mesh, the internal force is computed as follows: where B is the gradient matrix and N elem is the number of elements.
Since the central difference method is explicit, a finite stable time step size which is necessary for a sound wave to cross an element in the mesh must be below a critical value to provide numerical stability (Courant condition [32]): (15) where ∆x is the length of the smallest element in the mesh, c is the speed of sound in the material. For a solid material, the speed of sound is: where ρ is the material density, G is the shear modulus, and P (ρ, e ) is the equation of state. In Equation (17), the second term on the right hand side accounts for the stiffening effect due to the increase of internal energy as the material is compressed. For a fluid material, k = ρ 0 c 2 in which ρ 0 is the mass density and c is the sound velocity. The viscosity of the fluid material is generally ignored in the calculation of the speed of sound. For sloshing problems the pressure component of stress is much greater than the deviatoric part of the stress due to low the viscosity of the fluid, and the deviatoric stress is sometimes ignored.
Before the advection phase of operator split method, an interface-tracking algorithm is performed in order to compute accurately the material interfaces in the ALE elements containing several fluid materials.

STRESS EQUILIBRIUM AND INTERFACE TRACKING
After the Lagrangian phase is performed, either the stress tensor, pressure or deviatoric stress should be equilibrated, but most mixture theories equilibrate only pressure [30], the pressure equilibrium is a non-linear problem, which is complex and expensive to solve. Skipping the stress equilibrium phase is assuming an equal strain rate for both materials, which is incorrect. For most problems, the linear distribution based on volume fraction of the volumetric strain during the Lagrangian phase also leads to incorrect results. The volume distribution should be scaled by the bulk compression of the two materials in the element. For example, in an element containing air and water, the air, which is highly compressible, will absorb most of the volumetric strain. By assuming an equal strain rate or volumetric strain scaled on the volume fraction of the element, the water is forced to accept the same amount of strain as the air, and will undergo artificial high stresses.
The tracking of the material interfaces in the ALE elements containing several fluids can be performed by the VOF (Volume of Fluid) method or the Young method [35] which is attractive for solving a broad range of non-linear problems in fluid and solid mechanics such as, sloshing and explosion applications [19 and 36], because it allows arbitrary large deformations and enables free surfaces to evolve. In this method, different material occurrences are considered by their respective volume fractions on the element level. For multi-material elements, the volume fraction of one fluid satisfies: (18) The total stress by σ is weighed by volume fraction to get the fluid stress fields: (19) The Young or volume of fluid (VOF) method is originally developed to track an interface in elements containing two materials for two-dimensional problems. This method is adapted in this paper for the three-dimensional problems. In this method, the material layout is described solely by the volume fraction repartition of the fluid material in the ALE elements. Specifically, a straight line using the simple linear interface calculation (SLIC) technique of Woodward and Collela [37] approximates the interface in the cell described on Figure 2. Interfaces are initially drawn parallel to the element faces. Then nodal volume fraction is computed to each node based on the fraction volumes of elements that share the same node. This nodal volume fraction repartition determines the slope of the material interface inside the element. The normal vector to the interface inside the element is defined by: Multi-material eulerian cell where f is the nodal volume fraction. The position of the interface is then adjusted so that it divides the element into two volumes, which correctly matches the element volume fraction. The interface position is used to calculate the volume of the fluid flowing across cell sides.
As the X-advection, Y-advection and Z-advection are calculated in separate steps, it is sufficient to consider the flow across one side only.

ADVECTION PHASE
In the second phase, called advection or transport phase, the transportation of mass, momentum and energy across element boundaries are computed. This may be thought of as remapping the displaced mesh at the Lagrangian phase back to its initial position. The transport equations for the advection phase are: where is the difference between the fluid velocity , and the velocity of the computational domain , which will represent the mesh velocity in the finite element formulation. In some papers is referred as the convective velocity.
Equation (21) is solved successively for the conservative variables: mass, momentum and energy with initial condition φ 0 (x), which is the solution from the Lagrangian calculation of Equations (8)-(9) at the current time. In Equation (21), the time t is a fictitious time: in this paper, time step is not updated when solving for the transport equation. There are different ways of splitting the Navier-Stokes problems. In some split methods, each of the Stokes problem and transport equation are solved successively for half time step. The hyperbolic equation system (21) is solved for mass momentum and energy by using a finite volume method. Either a first order upwind method or second order Van Leer advection algorithm [38] can be used to solve Equation (21).

GOVERNING EQUATIONS FOR STRUCTURE Let
, the domain occupied by the structure, and let denote its boundary ( Figure 3). An updated Lagrangian finite element formulation is considered: the movement of the structure Ω s described by x i (t) (i = 1,2,3) can be expressed in terms of the reference coordinates X α (t) (α = 1,2,3) and time t, The momentum equation is given by Equation (24) in which is the Cauchy stress, ρ is the density, f is the force density, is acceleration and is the unit normal oriented outward at the boundary The solution of Equation (24)

FLUID-STRUCTURE INTERACTION
Fluid-structure interaction problems can be treated with either contact or coupling interface algorithms. Both algorithms are used to compute the interface forces applied from the fluid to the structure and conversely. These forces are imposed to the fluid and structure nodes in interface in order to prevent a node from passing through fluid-structure interface. For explicit methods, nodal forces at the fluid-structure interface are updated at each time step to account for interface forces. In contact problems, the slave and master meshes geometrically define the contact interface, whereas in the fluid-structure coupling method developed in this paper, the fluid coupling interface is defined by the material surface. From a mechanical point of view, the Euler-Lagrange coupling algorithm introduced in this paper is similar to penalty contact algorithm of Lagrangian analysis because the coupling method is mainly based on force equilibrium, and energy conservation. This method can be described as Eulerian contact. Therefore, in the following sections, a detailed description of the Euler-Lagrange coupling algorithm is presented together with the description of a regular penalty contact algorithm as well as the interface conditions for contact and coupling formulations.

INTERFACE CONDITIONS
Let us assume that a fluid Ω f and a structure Ω s given in Figure 4 are in contact. Therefore, the gap, d, normal to the common interface is zero. However, for clarity, the contacting surfaces are depicted in this figure separately. This sketch sums up Equations (2)-(7) for the fluid and the Equations (24)-(26) for the structure. In a fluid-structure interaction where and are normals at the contact point for Ω s and Ω f , respectively, and are the stress fields of the fluid and structure, respectively.
The impenetrability conditions for Ω f and Ω s can be stated as ∂ ∂ ∂ Figure 4 Interaction between Ω f and Ω s (the contacting surfaces are sketched separately for clarity).
It is convenient to express the condition Equation (28) in terms of which is the penetration rate (29) where and are the contact point velocities of the fluid and structure, respectively.
The condition Equation (29) expresses the fact that the fluid and structure must either remain in contact ( ) or separate ( ). It may appear inconsistent to speak of a penetration rate when impenetrability is an important condition on the solution. However, in many numerical methods, a small amount of interpenetration is allowed: d < 0. Then, the condition Equation (29) will not be observed with exactitude. In the penalty method, the impenetrability constraint is imposed as a penalty normal traction along the fluid-structure interface.

CONTACT ALGORITHMS FOR FLUID-STRUCTURE INTERACTION PROBLEMS
In contact algorithms, 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. For a fluid-structure interaction problem, the fluid nodes at the interface are considered as slave, and the structure elements as master. The structure is represented as Lagrangian whereas an ALE or Lagrangian mesh is used for the fluid.
There are two different approaches for solution of the contact problems. The first approach is the kinematic contact, constraining fluid and structure nodes to the same velocity. Kinematic contact conserves total momentum, but not total energy. The second approach, the penalty contact, is different from the previous one. The penalty method imposes a resisting force to the slave node, proportional to its penetration through the master element. This force is applied to both the slave node and the nodes of the master element in opposite directions to satisfy equilibrium. The resisting force applied to the slave nodes is computed according to following equation: The force applied to the nodes of the master element is scaled by the shape functions. Therefore, the force on one of the nodes of the master element: (31) where N i is the shape function at node i (i = 1, 2 in two dimensions, and i = 1, …, 4 in three dimensions) of the master element surrounding the slave node location, and d the penetration vector. In case the slave node coincides exactly with one of the master node ( Figure 5), node 1 for instance, we obtain: The coefficient k represents the stiffness of a spring. In fact, this method consists in placing springs between all penetrating nodes and the contact surface. The spring stiffness is given by Equation (33) in terms of the bulk modulus K of the master material, V the volume of the master element and A the area of the master segment: (33) where is scale factor for the interface stiffness, which satisfies .
Nevertheless, one of the problems encountered in contact applications is the high mesh distortion at the contact interface due to high fluid nodal displacements and velocities. Since the fluid nodes at the contact interface move in order to remain in contact with the Lagrangian structure, an ALE method is required to remesh the fluid domain. For small fluid mesh deformations, classical ALE methods described in [18] and [39] can be used but for large mesh distortion, ALE methods cannot be used for the mesh motion. In order to solve the problem, a rezoning or automatic remesh methods including the equipotential methods, simple and volume average, is required for the fluid domain, which is CPU time consuming. This method is based on interpolation algorithm, which is first-order accurate, non conservative and numerically dissipative ( [17] and [40]). An alternative algorithm to avoid fluid mesh distortion is to use an Euler-Lagrange coupling, an Eulerian formulation for the fluid and a Lagrangian formulation for the structure.

EULER-LAGRANGE COUPLING
Coupling methods are used to link the Lagrangian parts to the Eulerian. The coupling method described in this paper is based on the penalty method for contact algorithms. The fluid is represented by solving Navier-Stokes equations with an Eulerian formulation on a Cartesian grid that overlaps the structure, while the structure is discretised by a Lagrangian approach. Although the formulation can be extended to an ALE formulation, for simplicity, the numerical simulations have been restricted to an Eulerian formulation for the fluid. Since the mesh is fixed and the fluid material flows through the mesh, a larger Eulerian mesh than the physical fluid mesh is required. Slave node Fluid t = t n t = t n + ∆t Similar to the contact, there are two coupling methods available, the penalty-based method and the constraint-based method. The penalty-based coupling method is the preferred procedure, cause unlike the constraint-based coupling method, energy is conserved, even though there may be certain problems with stability.
In an explicit time integration problem, the main part of the procedure in the time step is the calculation of the nodal forces. After computation of fluid and structure nodal forces, coupling forces are added to the nodes that are on the fluid structure interface. For each structure node, a depth penetration is incrementally updated at each time step, using the relative velocity at the structure node, which is considered as a slave node, and the master node within the Eulerian element. The location of the master node is computed using the isoparametric coordinates of the fluid element. At time ( Figure 6) the penetration depth represents the amount the interface condition is violated, it is updated incrementally: The fluid velocity is the velocity at the master node location, interpolated from the nodes of the fluid element at the current time. For this coupling, the slave node is a structure mesh node, whereas the master node is not a fluid mesh node, it can be viewed as a fluid particle within a fluid element, with mass and velocity interpolated from the fluid element nodes using finite element shape functions. The vector represents the penetration depth of the structure inside the fluid during the time step, which is the amount the constraint is violated.
The coupling force acts only if penetration occurs, , where is built up by averaging normals of structure elements connected to the structure node. For clarity, the superscript of the penetration has been omitted, we will use d instead of , for the  penetration vector. 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 6.
Similarly to penalty contact algorithm, the coupling force is described by (35) where k represents the spring stiffness, and d the penetration. The force F in Equation (35) 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. At the structure coupling node, we applied a force (36) whereas for the fluid, the coupling force is distributed to the fluid element nodes based on the shape functions, at each node i (i = 1, …, 8), the fluid force is scaled by the shape function N i (37) Since the action-reaction principle is satisfied at the coupling interface and coupling force reduces fluid penetration into the structure. In contrast to the constrained method, where the interface condition is imposed and no interpenetration allowed, the penalty method allows some interpenetration. However, coupling forces proportional to the magnitude of the penetrations are applied to keep the penetration small. In the limit, defining a very high penalty stiffness, the penetrations approach zero, satisfying the interface condition. However, the coupling interface adds stiffness to the system affecting its eigenfrequency spectra. Hence, for numerical stability reasons, at a given time step size, there is an upper bound for the maximum possible penalty stiffness coefficient.
The main difficulty in the coupling problem comes from the evaluation of the stiffness coefficient k in Equation (35). 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. For fluid structure coupling, the spring stiffness is deduced from explicit contact with penalty method. This value is difficult to obtain unless numerical experiments are run systematically. For some industrial problems, automotive industry for instance, experimental tests are done prior to analysis to evaluate the stiffness value for penalty contact algorithm. In this paper, the stiffness k is based on stiffness used in explicit contact algorithms in [41]. k is given in terms 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 elements connected to the structure node: (38) However, to avoid numerical instabilities, a scalar factor is introduced as in the contact method. In order to avoid this anomaly, the force F in Equation (35) can be bounded by the contact force between two spheres, defined by Belytschko and Neal [42], which is given by Equation (39) (39) where M s is the mass of the slave or structure node, M m is the mass of the master fluid node interpolated from the fluid element nodes (Equation (40)), and is the relative velocity defined in Equation (34), is the current time step. The accuracy of the coupling also depends upon the number of coupling points used in the coupling procedure, besides the element nodes. Additional coupling points can be positioned on the structure, allowing a finer structure mesh only for the coupling. A compromise has to be made between to many additional coupling points, leading long computational times, or a few if no additional coupling points, resulting in the effect of artificial leakage through the structure.
As mentioned in the previous section, the coupling algorithm can be used for problems involving large mesh distortion that contact algorithm cannot handle since high mesh distortion decreases time step. However, problems related to fluid leakage through the structure may occur for high velocity problems. In such cases the fluid particle penetrates so deep within the structure that the coupling force in Equation (35), is not large enough to return it to the coupling interface. For general problems, one solution to this problem is to reduce the time step; but for some problems involving highly compressible gas, the expression Equation (35) needs to be modified: the spring system is nonlinear and the stiffness k should be penetration dependent: (41) However, fluid leakage through the structure is a very difficult problem to solve. Most of the research in fluid structure interaction has focused on developing numerical algorithms that prevent leakage and conserve energy.

SLOSHING IN A RIGID TANK
Sloshing in liquid containment tanks has caused serious damage on tank wall and roof during many earthquakes. The correct description of these forces is very important to provide structural safety of tanks. In order to evaluate sloshing effects in tanks, a 3D rectangular rigid tank model subjected to harmonic motion is analyzed using the coupling algorithm presented in this paper. Time history of the free surface wave height and pressure obtained at specific locations are compared with the existing experimental, analytical and numerical results. The experimental study carried out by Liu and Lin [18] is regarded as reference solution to evaluate the numerical findings. The tank dimensions and material properties are determined in accordance with this reference model. The analytical method developed by Faltinsen [9] which is based on the linear potential flow theory is used to compare the nonlinear effects of fluid sloshing with those of linear. In his work, tank problem with sloshing liquid is treated as an initial boundary-value problem. Liquid is assumed homogeneous, inviscid, irrotational and incompressible. The wave amplitudes are very small in comparison with the wavelengths and depths. These assumptions permit representation of fluid flow inside the tank with Laplace equation. Furthermore, the algorithm based on the ALE method detailed by Aquelet et. al [19] is applied to the same tank problem and results are evaluated. In this method, the fluid structure interface is described by the mesh nodes which are sheared by the fluid and the structure at the interface. The fluid region is treated on a moving mesh using an ALE formulation whereas the structure is discretized with a deformable mesh using a Lagrangian formulation. This method is suitable for fluid structure interaction problems where structure deforms considerably moderate. For these kind of problems, ALE mesh can accord the deformations of Lagrangian mesh at the fluid structure interface and stable time step size is high enough for explicit time integration algorithms to continue.
The three-dimensional rigid rectangular tank used in the numerical model has a width of 0.57 m, breath of 0.31 m and total height of 0.30 m and is filled with water (ρ = 1000 kg/m 3 ) up to a height of 0.15 m. The interior liquid is discretized with uniform mesh (Figure 7). Hydrostatic pressure field is generated increasing gradually until 1 sec. The model is fixed the base and two harmonic motions with resonant and non-resonant frequencies are applied to the model. The first sloshing mode circular frequency is obtained as ω o = 6.0578 rad/s by the following equation: (42) where, n, L and h represent mode number, tank length and fluid depth, respectively. Therefore, the excitation frequencies are determined as ω = 0.583 ω o and ω o for nonresonance and resonance cases, respectively. The displacement amplitudes of the horizontal harmonic excitations are 0.005 m for both cases. The time history response of free surface elevation is measured at three locations which are near left (i.e. x = −0.265) and right (i.e. x = 0.265) ends of tank and at the middle of the free surface (i.e. x = 0).
For non-resonant frequency motion, the numerical solution of sloshing by the proposed method is in a quite acceptable agreement with the reference solution and analytical formulation in terms of elevation of free surface. Figure 8  For negative (downward) wave amplitudes, numerical results are more consistent than those of analytical, since numerical method takes into account nonlinear sloshing behavior. As it is expected, the wave height is almost zero at the middle of the free surface for each method.
For the resonant frequency case, the wave height increases continuously over time for all solution types at the near left and right end of the tank. The comparison of three solution methods shows that analytical study overestimates negative surface amplitudes, whereas it underestimates the positive ones ( Figure 9). Numerical and experimental results are highly consistent in terms of peak level timing, shape and amplitude of sloshing wave. The free surface displacement time histories obtained from numerical and experimental studies show that the positive (upward) sloshing wave amplitudes are always larger than the negative (downward) ones. This phenomenon is a classical indication of a nonlinear behavior of sloshing and caused by suppression effect of the tank base on waves with negative amplitude. Although the gravity effects exist for upward and downward fluid motion, the downward motion of fluid is blocked by the tank bottom. The ratio of positive amplitude to absolute negative amplitude increases as the fluid depth decreases. This phenomenon can not be observed from analytical solution because it is derived under linearized assumptions.
There is a strong correlation between two numerical methods where the fluid is represented by solving Navier-Stokes equations with an ALE formulation and structure is treated on a deformable mesh using a Lagrangian formulation. The only difference of two methods is the connectivity of both material nodes at the fluid structure interface. In the coupling method, fluid and structure have independent nodes whereas in the numerical method proposed by Aquelet et al. [19] fluid and structure have common nodes at the interface of the two materials. For both resonant and non-resonant frequency cases, free surface time history results for both numerical methods are almost the same in terms of peak level timing, shape and amplitude of sloshing wave (Figures 10 and 11).
Since the pressures were not measured during the experimental study, the analytical results are considered as reference solution. Yet, it should be bore in mind that the analytical method does not include the nonlinear effects. Pressure (including hydrostatic pressure) results are observed at three specific locations which are at the two edge of the tank wall, 0.01 m above the base, and at the middle of the tank base plate. There is a high frequency oscillation region at the beginning of the pressure plots, but oscillations disappear after 4 seconds. For both the resonance and non-resonance cases, the presented numerical method predicts the pressure higher than the analytical and other numerical method at the bottom of the tank wall (x = -0.265 and 0.265 m). However, the pressure measured at the middle of the tank base (x = 0.0 m) is exactly the same for both numerical methods. At this location, the analytical method gives a single value throughout the analysis (Figures 12 and 13). For resonant frequency loading case, pressures obtained from analytical study at the right and left walls of the tank increase continuously over time. However, pressure observed from numerical models fluctuates between the same negative and positive values after 8 seconds. This verifies that analytical method is not reliable for resonant frequencies where nonlinear sloshing behavior is extremely dominant.

APPLICATION TO A CYLINDRICAL FLEXIBLE LIQUID STORAGE TANK
The applicability of coupling algorithm presented in this paper for seismic tank problems is verified with the results of the existing experimental study carried out by Manos and Clough [26]. The analyses are performed for the same tank dimensions and material properties which are used in the experimental study.
The tank is made of aluminum with a density of 2700 kg/m 3   placed on the top of the second tank shell course which has a thickness of 0.0013 m. Arbitrary Lagrangian Eulerian (ALE) description of the liquid-structure interface is employed in order to enforce compatibility between structure and liquid elements. In the numerical simulations, both material and geometric nonlinearities are considered in order to accurately determine stress, strain and strain rate distributions throughout the tank and fluid. Since the free surface motion of fluid is nonlinear in reality, single plane of geometric and loading symmetry in the structure are not used in the analyses; the response of the whole system is determined. The resulting discretized model which has 15162 nodes, 1120 shell elements and 13056 solid elements is given in Figure 14. The tank base is constrained in five degrees-of-freedom using single point constraints on all nodes. A vertical acceleration field of a 1 g is applied to give the correct hydrostatic pressure in the fluid. Nonlinear dynamic time history analyses are performed under unidirectional horizontal earthquake motion recorded during the 1940 El Centro earthquake with 0.50 g peak acceleration, which is scaled with regard to time by 1/ (Figure 15). Numerical results obtained by the presented method are highly consistent with the experimental results in terms of peak level timing, shape and amplitude of sloshing wave. 14 16 18 20 x = 0.265; y = 0.01 Figure 12 Comparisons of the time histories of pressure for the present numerical method, numerical method proposed by Aquelet et al. [19] and analytical data (non-resonant case).
on the loading axis. The maximum sloshing amplitudes at x = -1.72 m are measured as 0.08 m for both downward and upward directions from the experimental study whereas, it is obtained as 0.08 and 0.045 m from the numerical solution for positive and negative amplitudes, respectively. At the second measurement location, numerical method gives 0.07 m for both positive and negative of sloshing amplitudes while negative and positive sloshing amplitudes from experimental study are 0.1 and 0.05 m, respectively. It can be noticed from these figures that the numerical and experimental models lead to a relatively accurate description of the water sloshing in terms of the displacement of free surface for the input earthquake motion. Therefore, it can be justified that the presented method is reliable for the analysis of sloshing inside a cylindrical tank subjected to earthquake excitation. For this earthquake motion, as similar to experimental findings, numerical simulation verifies the existence of the high frequency, resonant out-of-round distortions of the circular cross-section of the tank in response to earthquake motion. Since fluid and structure meshes are independent these distortions do not cause any large mesh deformation problems. Plastic deformations do not develop during both experimental and numerical studies and tank behaves in an elastic manner.

CONCLUSIONS
In this study, a fully non-linear numerical coupling algorithm based on Arbitrary Lagrangian Eulerian (ALE) method is developed to solve the three-dimensional fluid-structure interaction problems. In this method, Navier Stokes equations are employed to represent the motion of the fluid domain which is treated on a fixed or moving mesh using an ALE formulation. An independent deformable mesh using a Lagrangian formulation is employed for the structure. The volume-of-fluid (VOF) method is adopted to track the free surface motion. The momentum advection algorithm is applied to transfer the materials across element boundaries. A second order accurate time marching procedure based on central difference method is used to advance the solution in time. In order to validate the applicability of the numerical algorithm to sloshing tank problems, simulations on a rigid rectangular tank model are carried out for non-resonant and resonant loading cases. The results of the numerical algorithm are compared with those of analytical method established on potential flow theory and experimental data in terms of sloshing wave elevation. It is obviously verified that the free surface profiles obtained from experimental and numerical studies perfectly match each other for both resonance and nonresonance loading conditions, although analytical results highly deviate from experimental results for loading with resonant frequency. The pressures are obtained at specific locations and compared with linear analytical results. Though pressures obtained by the analytical method continuously amplify over the time for resonance loading case, numerical method gives bounded results after 8 seconds. In the sequential numerical analysis, the sloshing behavior of a flexible cylindrical tank observed during the laboratory tests is correlated with numerical simulations in terms of sloshing wave height. Comparison results showed that the present numerical algorithm is reliable and useful for sloshing problems in practice for every frequency range of external excitation.