Int. Jnl.

of Multiphysics Volume 3 · Number 3 · 2009 307

ALE and Fluid Structure Interaction

for Sloshing Analysis
Z. Ozdemir1,4, M. Moatamedi2, Y.M. Fahjan3 and M. Souli4
1BogaziciUniversity, Kandilli Observatory, Cengelkoy, Istanbul, TURKEY
2Cranfield University, Cranfield, Bedfordshire, MK43 0AL UK
3Gebze Institute of Technology (GYTE), Cayirova Campus,

Kocaeli, TURKEY
4Université de Lille, Laboratoire de Mécanique de Lille, CNRS 1441,

Bd Paul Langevin Villeneuve d’Ascq, FRANCE

ABSTRACT
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 fluid-
structure 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.

Keywords: Finite element method, fluid-structure interaction, earthquake,

liquid containment tanks, penalty based coupling method
308 ALE and Fluid Structure Interaction for Sloshing Analysis

1. 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
methods regarded the tank as rigid ([1], [2], [3] and [4]), earthquake experiences (the 1964
Alaska, the 1971 San Fernando and the 1979 Imperial County California earthquakes)
verified that the fluid-structure interaction effects between two materials cause amplification
of the peak ground acceleration several times. Furthermore, the impact of large amplitude
nonlinear sloshing waves on tank wall and roof generates additional hydrodynamic forces
which can affect the overall response of the tank. Therefore, the analysis methodology used
for the solution of tank problems has to be capable of handling accurately the highly
nonlinear behavior of tank itself, contained liquid and coupling effects of fluid and structure.
The first approximate mathematical solution of impulsive hydrodynamic effects on a fully
anchored rigid cylindrical tank under horizontal motion was provided by Jacobsen [1]. In this
study, the irrotational motion of invicid and incompressible fluid inside a rigid cylindrical
container was represented with Laplace equation which satisfies specified boundary
conditions. With the same assumptions, Housner [2, 3 and 4] split hydrodynamic effects due
to seismic motions in two parts namely impulsive effects caused by the liquid portion which
accelerates with the tank and convective effects caused by liquid part which participates
sloshing motion at the free surface. He represented the impulsive and convective masses on a
simple spring mass model which constructs the basis of the analysis method of the current
design codes. Veletsos and Yang [5] added the effect of the convective liquid motion to the
Jacobsen’s approach and expressed the solution of the Laplace’s equation as the sum of the
impulsive and convective components. In order to take into account flexibility of a rigidly
anchored tank, Veletsos and Yang [5] applied Flügge’s shell theory with a Rayleigh-Ritz type
procedure using the natural modes of vibration of uniform cantilever beams. Haroun and
Housner [6 and 7] used boundary integral theory to model fluid region and ring shaped finite
elements for tank shell to develop reliable method for analyzing the dynamic behavior of
deformable cylindrical tanks. Veletsos [8] improved Housner’s mechanical analogue in order
to take into account the effect of tank wall flexibility. Faltinsen [9] derived a linear analytical
solution for liquid sloshing in a horizontally excited 2-D rectangular tank considering
damping due to viscous effects. Fischera and Rammerstorferb [10] investigated analytically
the overall effect of pressure generated by interaction forces between sloshing and the wall
motion modifying the free surface boundary conditions. In the recent years, it was realized
that the behavior of unanchored tanks is quite different than that of rigidly anchored. Malhotra
and Veletsos [11] studied extensively the uplifting behavior of the bottom plate of an
unanchored tank by idealizing the base plate as a uniformly loaded semi-infinite, prismatic
beam that rests on a rigid foundation. The effect of elastic end constraints, the influence of the
axial forces associated with large deflections and the effect of plastic yielding in the beam
were considered. A sequel to this work, the effects of uplifting of bottom plate of the tank on
the entire tank-liquid system was thoroughly investigated ( [12] , [13] and [14] ).
The numerical techniques, especially the finite element method, have been widely used
for replicating actual physical behavior of not only the tank itself but also the contained fluid
with greater accuracy. El-Zeiny [15] developed a finite element program, which uses an
updated Eulerian-Lagrangian description of the liquid-structure interface, in order to analyze
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 309

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.
310 ALE and Fluid Structure Interaction for Sloshing Analysis

2. 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.

2.1. 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
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 311

h (t)
∂Ωf1

Ωf g (t)

∂Ωf2

Figure 1 Fluid domain and boundary condition description.

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.

2.2. 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:
r r
d f ( X , t ) ∂ f ( x , t ) r r uuuuur r
= + ( v − w). grad f ( x , t ) (1)
dt ∂t

r r r
where X is the Lagrangian coordinate, x is the ALE coordinate, v is the particle velocity
r
and w 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 Ω f ∈ R3 represent the domain occupied by the fluid particles, and let ∂ Ω f 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:

∂ρ r r r
+ ( v − w) grad ( ρ ) + ρ div ( v ) = 0 (2)
∂t

r
∂v r r r uuur r
ρ + ρ ( v − w) . grad ( v ) = div (σ ) + f (3)
∂t
312 ALE and Fluid Structure Interaction for Sloshing Analysis

∂e r r uuuuur r rr
ρ + ρ ( v − w). grad ( e) = σ : grad ( v ) + f .v (4)
∂t

where ρ is the density and σ is the total Cauchy stress given by:

r r
σ = − p.Id + µ ( grad (v ) + grad ( v )T ) (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 ∂ Ω1f . The inflow
boundary condition is:

r r
v = g (t ) on ∂ Ω1f (6)

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.

r r
σ . n = h (t ) on ∂ Ω 2f (7)

r
where, n is the outward unit normal vector on the boundary.
One of the major difficulties in time integration of the ALE Navier-Stokes Equations
r r
(2)–(4) is due to the nonlinear term related to the relative velocity ( v − w ) 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,
r since,
r the problem is formulated in the spatial coordinate, the reference frame
is fixed: w = 0 . 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
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 313

velocity, pressure and internal energy due to external and internal forces are computed. The
equilibrium equations for the Lagrangian phase are:
r
dv uuur ur
ρ = div (σ ) + f (8)
dt

de r ur r
ρ = σ : grad ( v ) + f .v (9)
dt

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]:

ρ J = ρ0 (10)

where ρ0 is the initial density and J is the volumetric strain given by the Jacobian:

 ∂x 
J = det  i  (11)
 ∂ X j 

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:

n+ 12
x n+1 = x n + ∆ t n u (12)

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:

n+ 12 n− 12
u =u + ∆t n . M −1.( Fext
n
+ Fint
n
) (13)

where Fint is the internal force vector and Fext 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:
314 ALE and Fluid Structure Interaction for Sloshing Analysis

Nelem
Fint = ∑ ∫ Bt .σ .dv (14)
k =1 k

where B is the gradient matrix and Nelem 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]):

∆x
∆t ≤ (15)
c

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:

4
G+k
(16)
c 2
= 3
ρ0

∂P P ∂P
k = ρ0 + (17)
∂ρ ρ ∂e

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 = ρ0c 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.

2.3. 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.
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 315

gradf
Fluid 2 n=
gradf

Material
interface

Fluid 1

Multi-material eulerian cell

Figure 2 Interface tracking in a fluid cell with volume of fluid (VOF) method.

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:

Vf ≤1 (18)

The total stress by σ is weighed by volume fraction to get the fluid stress fields:

σ f = σ .V f (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:
uuuuur
r grad f
n = uuuuur (20)
grad f
316 ALE and Fluid Structure Interaction for Sloshing Analysis

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.

2.4. 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:

∂φ r uuuuur (21)
+ c.grad (φ ) = 0
∂t
r
φ ( x , 0) = φ0 ( x ) (22)

r r r r
where c = v − w is the rdifference between the fluid velocity v , and the velocity of the
computational domain w , whichr will represent the mesh velocity in the finite element
formulation. In some papers c 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).

2.5. GOVERNING EQUATIONS FOR STRUCTURE

Let Ω S ∈ R3 , the domain occupied by the structure, and let ∂ Ω denote its boundary
S

(Figure 3). An updated Lagrangian finite element formulation is considered: the movement
of the structure Ω s described by xi(t) (i = 1,2,3) can be expressed in terms of the reference
coordinates Xα (t) (α = 1,2,3) and time t,

x i = x i ( Xα , t ) (23)

The momentum equation is given by Equation (24) in which σ is the Cauchy stress,
r r
ρ is the density, f is the force density, dv is acceleration and n is the unit normal
dt
oriented outward at the boundary ∂ Ω S
r
dv uuur ur
(24)
ρ = div (σ ) + f
dt
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 317

τ (t )
∂Ωs1

Ωs D (t )

∂Ωs2

Figure 3 Structure domain and boundary condition description.

The solution of Equation (24) satisfies the displacement boundary condition Equation
(25) on the boundary ∂ Ω1S and the traction boundary condition Equation (26) on the

boundary ∂ Ω 2
S

r r r
x ( X , t ) = D(t ) on ∂ Ω1S (25)

r r
σ . n = τ (t ) on ∂ Ω 2S (26)

3. 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.

3.1. 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
318 ALE and Fluid Structure Interaction for Sloshing Analysis

τ (t)
∂ Ωs1

Ωs D (t)

s
∂ Ω2

n2
v2

vf
h (t)
nf
∂Ωf1

Ωf g (t)

∂ Ωf2

Figure 4 Interaction between Ω f and Ω s (the contacting surfaces are sketched

separately for clarity).

problem, two conditions applied on the interface ∂ Ω C2 common to ∂ Ω 2S and ∂ Ω 2f

are added to the previous equations: the fluid–structure traction conditions and the
impenetrability condition. The traction condition observes the balance of momentum across
the fluid–structure interface. Since this interface has no mass, the sum of traction forces on
the fluid and structure must vanish (Newton’s third law). Since a frictionless model is
considered, the tangential tractions vanish too. Thus the traction condition is simply the
normal traction relation:

S r f r
σ . n S + σ . n f = 0 on ∂ Ω C2 (27)
rS rf S f
where n and n 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

Ω f ∩ Ω S = 0 on ∂ Ω C2 (28)
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 319

& which is the

It is convenient to express the condition Equation (28) in terms of d,
penetration rate

r r r r r r r
d& = v S .n S + v f .n f = ( v S − v f ).n S ≤ 0 on ∂Ω C2 (29)

rS rf
where v and v 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 ( d& = 0 ) or separate ( d& < 0 ). It may appear inconsistent to speak of a
penetration rate d, & 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.

3.2. 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:

FS = − k ⋅ d (30)

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:

Fmi = N i ⋅ k ⋅ d (31)

where Ni 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:

Fm1 = k ⋅ d and Fm2 = Fm3 = Fm4 = 0 (32)

320 ALE and Fluid Structure Interaction for Sloshing Analysis

Slave node
Master node
Structure Ff = −k ⋅ d

d k

Fluid Master node Slave node

Ff = k ⋅ d

t = tn t = tn + ∆t

Figure 5 Sketch of the contact algorithm.

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:

K A2 (33)
k = pf ⋅
V

where p f is scale factor for the interface stiffness, which satisfies 0 ≤ p f ≤ 1 .

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.

3.3. 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.
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 321

Void Void
Structure
Ff = −k ⋅ d

Slave
node
Fluid
Fluid
element Fluid
Fs = k ⋅ d
Fluid
node
Master node
(fluid particle) d
t= tn t = tn + ∆t

Figure 6 Sketch of the coupling algorithm.

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
r
each structure node, a depth penetration d is incrementally updated at each time step, using
r r
the relative velocity ( vS − v f ) 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 t = t n (Figure 6)
r
d n the penetration depth represents the amount the interface condition is violated, it is
updated incrementally:

r r r n+ 1 r n+ 1
d n +1 = d n + ( vS 2 − v f 2 ) ∆ t (34)

r
The fluid velocity v f 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
r
within a fluid element, with mass and velocity interpolated from the fluid element nodes
using finite element shape functions. The vector d n represents the penetration depth of the
structure inside the fluid during the time step, which is the amount the constraint is violated.
r r r
The coupling force acts only if penetration occurs, ns d n < 0 , where ns is built up by
averaging normals of structure elements connected to the structure node. For clarity, the
rn
superscript of the penetration has been omitted, we will use d instead of d , for the
322 ALE and Fluid Structure Interaction for Sloshing Analysis

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

F = k⋅d (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

FS = − F (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 Ni

F fi = N i ⋅ F (37)

8
Since ∑ F fi = F , the action–reaction principle is satisfied at the coupling interface
i=1
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
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 323

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:

K A2 (38)
k = pf ⋅
V

However, to avoid numerical instabilities, a scalar factor p f , 0 ≤ p f ≤ 1, 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)
r r
M S M m ( vS − v f )
F≤ (39)
∆ t ( MS + Mm )

8
Mm = ∑ N i ⋅ M mi (40)
i =1

where Ms is the mass of the slave or structure node, Mm is the mass of the master fluid
r r
node interpolated from the fluid element nodes (Equation (40)), and ( vS − v f ) is the relative
velocity defined in Equation (34), ∆ t 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:

F = k (d ) ⋅ d (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.
324 ALE and Fluid Structure Interaction for Sloshing Analysis

4. 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/m3)
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:

( 2 n + 1) π  ( 2 n + 1) π 
ω n2 = g tanh  h (42)
L  L 

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 non-
resonance 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 presents the time history response
of free surface elevations at three measurement location (i.e. x = −0.265, 0 and 0.265 m) that
is extended to 20 s. The maximum sloshing amplitude at x = +0.265 m is measured as 0.0069 m
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 325

Time = 6.68 Tank Fluid mesh

Y
X

Figure 7 Finite element model of the rigid rectangular tank with fluid.

from the numerical solution whereas, it is obtained as 0.0067 m and 0.0064 m from the
analytical formulation and experimental data, respectively. The minimum free surface
elevation is obtained as 0.0065 m, 0.0066 m and 0.0066 m from experimental, analytical and
numerical methods, respectively. In both cases, the numerical results are quite close to the
experimental ones. 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
326 ALE and Fluid Structure Interaction for Sloshing Analysis

Surface displacement (m) 0.01

x=0 ALE coupling
Analytical solution
0.005 Experimental

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

0.01
Surface displacement (m)

x = −0.265
0.005

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

0.01
Surface displacement (m)

x = 0.265
0.005

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

Figure 8 Comparisons of the time histories of free surface elevation for the
present numerical method, the analytical solution and experimental data (non-
resonant case).

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
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 327

0.15

Surface displacement (m)

x=0 ALE coupling
0.1 Analytical solution
Experimental
0.05
0
−0.05
−0.01
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

0.15
Surface displacement (m)

x = −0.265
0.1
0.05
0
−0.05
−0.01
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

0.15
Surface displacement (m)

x = 0.265
0.1
0.05
0
−0.05
−0.01
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

Figure 9 Comparisons of the time histories of surface elevation for the present
numerical method, the analytical solution and experimental data (resonant case).

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.
328 ALE and Fluid Structure Interaction for Sloshing Analysis

Surface displacement (m) 0.01

x=0 ALE coupling
ALE common nodes
0.005

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

0.01
Surface displacement (m)

x = −0.265 ALE coupling

ALE common nodes
0.005

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

0.01
Surface displacement (m)

x = 0.265 ALE coupling

ALE common nodes
0.005

−0.005

−0.01
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

Figure 10 Comparisons of the time histories of surface elevation for the present
numerical method and numerical method proposed by Aquelet et al. [19]
(non-resonant case).

5. 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/m3, elastic modulus of 71.0 GPa
and yield stress of 100 MPa. Model has a radius of 1.83 m and a total height of 1.83 m. Water
(ρ = 1000 kg/m3) is filled up to a height of 1.53 m. The thicknesses of the aluminum bottom
plate and the shell section nearest to the bottom is 0.002 m. L shaped steel wind girder is
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 329

0.15
Surface displacement (m)
x=0 ALE coupling
0.1 ALE common nodes
0.05
0
−0.05
−0.1
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

0.15
Surface displacement (m)

x = −0.265
0.1
0.05
0
−0.05
−0.1
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

0.15
Surface displacement (m)

x = 0.265
0.1
0.05
0
−0.05
−0.1
−0.15
0 1 2 3 4 5 6 7
Time, t(sec)

Figure 11 Comparisons of the time histories of surface elevation for the present
numerical method and numerical method proposed by Aquelet et al. [19]
(resonant case).

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.
330 ALE and Fluid Structure Interaction for Sloshing Analysis

2500
x = 0.0; y = 0.0 ALE coupling
2000 ALE common nodes
Pressure (N/m2)

Analytical
1500

1000

500

0
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

3000
x = −0.265; y = 0.01 ALE coupling
2500 ALE common nodes
Pressure (N/m2)

Analytical
2000
1500
1000
500
0
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

3000
x = 0.265; y = 0.01 ALE coupling
2500 ALE common nodes
Pressure (N/m2)

Analytical
2000
1500
1000
500
0
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

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).

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 uni-
directional 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/ 3 (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.
Figures 16 presents the time history response of free surface elevations at two measurement
locations (i.e. x = –1.72 and 1.72 m) that is extended to 7 sec. These locations are situated
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 331

2500
x = 0.0; y = 0.0 ALE Coupling
coupling
2000 ALE Common
common nodes
Nodes
Pressure (N/m2)
Analytical
1500
1000

500

0
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

3500
3000 x = −0.265; y = 0.01
Pressure (N/m2)

2500
2000
1500
1000
500
ALE coupling
Coupling
0 ALE common
Commonnodes
Nodes
−500 Analytical
−1000
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

3500
3000 x = 0.265; y = 0.01
Pressure (N/m2)

2500
2000
1500
1000
500
ALE coupling
Coupling
0 ALE common
Commonnodes
Nodes
−500 Analytical
−1000
0 2 4 6 8 10 12 14 16 18 20
Time, t(sec)

Figure 13 Comparisons of the time histories of pressure for the present numerical
method, numerical method proposed by Aquelet et al. [19] and analytical data
(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.
332 ALE and Fluid Structure Interaction for Sloshing Analysis

Time = 7.09 Fluid mesh

z Tank
y x

Figure 14 Finite element model of the cylindrical tank.

0.6

0.4
Acceleration, g

0.2
0
−0.2
−0.4
−0.6
0 1 2 3 4 5 6 7
Time, sec

Figure 15 Input motion used in the experimental study of cylindrical tank

(Manos and Clough [26]).

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.
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 333

0.1
ALE coupling
0.08 Experimental (EERC-82/07)

Free surface displacement (m) 0.06

0.04

0.02

−0.02

−0.04

−0.06

−0.08

−0.1
0 1 2 3 4 5 6 7
Time (sec)

0.1
ALE coupling
0.08 Experimental (EERC-82/07)

0.06
Free surface displacement (m)

0.04

0.02

−0.02

−0.04

−0.06

−0.08

−0.1
0 1 2 3 4 5 6 7
Time (sec)

Figure 16 Comparisons of the time histories of surface elevation for the present
numerical method and experimental data (left: x = –1.72 m; right: x = 1.72 m).
334 ALE and Fluid Structure Interaction for Sloshing Analysis

6. 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 non-
resonance 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.

ACKNOWLEDGEMENT
This work is supported by Engineering Research Group of the Scientific and Technological
Research Council of Turkey (TUBITAK-MAG) through contract numbers 108M607.

REFERENCES
[1] L.S. Jacobsen, 1949, Impulsive hydrodynamics of fluid inside a cylindrical tank and of fluid
surrounding a cylindrical pier. Bulletin of the Seismological Society of America, 39, 3, 189–204, July.
[2] G.W. Housner, 1954, Earthquake pressures on fluid containers. 8th Technical Report under Office
of Naval Research, California Institute of Technology, Pasadena, California, August.
[3] G.W. Housner, Dynamic pressures on accelerated fluid containers. Bulletin of the Seismological
Society of America, 47, 1, 15–35, January, 1957.
[4] G.W. Housner, The dynamic behavior of water tanks. Bulletin of the Seismological Society of
America, 53, 2, 381–387, February, 1963.
[5] A.S. Veletsos, J.Y. Yang, Earthquake response of liquid storage tanks. Advances in Civil
Engineering Through Engineering Mechanics, Proceedings of the Engineering Mechanics
Division Specialty Conferences, ASCE, Raleigh, North Carolina, 1–24, 1977.
[6] Haroun M. A. and Housner, G. W., 1981a, “Seismic Design of Liquid Storage Tanks”, Journal of
the Technical Councils of ASCE, Vol. 107, No. TC1, April.
[7] Haroun MA and Housner GW, 1981b, “Earthquake response of deformable liquid storage tanks”.
Journal of Applied Mechanics, ASME, Vol. 48, pp. 411–418., June.
Int. Jnl. of Multiphysics Volume 3 · Number 3 · 2009 335

[8] Veletsos A.S., 1984, “Seismic Response and Design of Liquid Storage Tanks, in Guidelines for
the Seismic Design of Oil and Gas Pipeline Systems”, ASCE, p. 255–370, 443–460.
[9] Faltinsen, O.M., “A numerical nonlinear method of sloshing in tanks with two-dimensional flow”,
J. Ship Res. 22 (1978) 193–202.
[10] F.D. Fischera, and F.G. Rammerstorferb, 1999, “A Refined Analysis of Sloshing Effects in
Seismically Excited Tanks”, International Journal of Pressure Vessels and Piping, Vol. 76,
pp. 693–709.
[11] Malhotra, P. K. and A. S. Veletsos, 1994a, “Beam Model for Base-Uplifting Analysis of
Cylindrical Tanks”, Journal of Structural Engineering, ASCE, Vol. 120, No. 12, pp. 3471-3488,
December.
[12] Malhotra, P. K. and A. S. Veletsos, 1994b, “Uplifting Analysis of Base Plates in Cylindrical
Tanks”, Journal of Structural Engineering, ASCE, Vol. 120, No. 12, pp. 3489–3505, December.
[13] Malhotra, P. K. and A. S. Veletsos, 1994c, “Uplifting Response of Unanchored Liquid-Storage
Tanks”, Journal of Structural Engineering, ASCE, Vol. 120, No. 12, pp. 3525–3547, December.
[14] Malhotra, P. K., 1995, “Base Uplifting Analysis of Flexibly Supported Liquid-Storage Tanks”,
Earthquake Engineering and Structural Dynamics, Vol. 24, Issue 12, pp. 1591–1607.
[15] El-Zeiny, A., 1995, “Nonlinear Time-Dependent Seismic Response of Unanchored Liquid
Storage Tanks”, PhD Dissertation, Department of Civil and Environmental Engineering,
University of California, Irvine.
[16] Barton DC and Parker JV., 1987, “Finite element analysis of the seismic response of anchored and
unanchored liquid storage tanks”, Earthquake Engineering and Structural
Dynamics;15(3):299–322.
[17] Souli, M., Ouahsine, A. and Lewin, L., “ALE formulation for fluid-structure interaction
problems”, Computer Methods in Applied Mech. and Eng. 190 (2000) 659–675.
[18] Souli, M and Zolesio, J P. Arbitrary Lagrangian-Eulerian and free surface methods in fluids
mechanics, Computational Methods In Applied Mechanics and Engineering, 2001 191 451–466.
[19] Aquelet, N., Souli, M., and Olovson L., “Euler Lagrange coupling with damping effects:
Application to slamming problems”, Computer Methods in Applied Mechanics and Engineering,
Vol. 195 , pp 110–132, (2005).
[20] Chen, Y. H., Hwang, W. S. and Ko, C. H., 2007, “Sloshing Behaviours of Rectangular and
Cylindrical Liquid Tanks Subjected to Harmonic and Seismic Excitations”, Earthquake
Engineering and Structural Dynamics, Vol. 36, pp.1701–1717.
[21] Liu, D. and Lin, P., “A numerical study of three-dimensional liquid sloshing in tanks”, Journal of
Computational Physics, 227 (2008), 3921–3939.
[22] Mitra, S., Upadhyay, P. P., and Sinhamahapatra, K. P., 2008, “Slosh Dynamics of Inviscid Fluids
in Two-Dimensional Tanks of Various Geometry Using Finite Element Method”, International
Journal for Numerical Methods in Fluids, Vol. 56, pp.1625–1651.
[23] M. S. Razzaghi and S. Eshghi, 2004, “Behavior of Steel Oil Tanks due to Near-Fault Ground
Motion”, 13th World Conference on Earthquake Engineering, Vancouver, B.C., Canada, August 1–6.
[24] Kana, D. D., 1979, “Seismic Response of Flexible Cylindrical Liquid Storage Tanks”, Nuclear
Engineering and Design, Vol. 52, pp. 185–199.
[25] Clough, D. P., 1977, Experimental evaluation of seismic design methods for broad cylindrical
tanks, UCB/EERC-77/10, PB-272 280.
336 ALE and Fluid Structure Interaction for Sloshing Analysis

[26] Manos and Clough, 1982, “Further study of the earthquake response of a broad cylindrical liquid-
storage tank model” UCB/EERC-82/07, Earthquake Engineering Research Center, University of
California, Berkeley, EERC-82/07.
[27] Manos, G. C., 1986, “Dynamic response of a broad storage tank model under a variety of
simulated earthquake motions.” Proc. ~ 3rd U.S. Nat. Conf. on Earthquake Engrg., Earthquake
Engineering Research Institute, E1 Cerrito, Calif., 2131–2142.
[28] Tanaka, Motoaki, Sakurai, Ishida, Tazuke, Akiyama, Kobayashi and Chiba, 2000, “Proving Test
of Analysis Method on Nonlinear Response of Cylindrical Storage Tank Under Severe
Earthquakes”, Proceedings of 12th World Conference on Earthquake Engineering (12 WCEE),
Auckland, New Zealand.
[29] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John
Wiley & Sons, New York, 2000.
[30] D.J. Benson, A mixture theory for contact in multi-material eulerian formulations, Comput. Meth.
Appl. Mech. Eng. 140 (1997) 59–86.
[31] T.J.R. Hughes, W.K. Liu and T.K. Zimmerman, “Lagrangian Eulerian finite element formulation
for viscous flows”, J. Computer Methods Applied Mechanics and Engineering, 29 (1981)
329–349.
[32] Benson, D.J., “Computational Methods in Lagrangian and Eulerian Hydrocodes”, Computer
Methods in Applied Mech. and Eng. 99 (1992) 235-394.
[33] Flanagan, D.P. and Belytschko, T., “A Uniform Strain Hexahedron and Quadrilateral and
Orthogonal Hourglass Control”, Int. J. Numer. Meths, Eng. 17,679–706 (1981).
[34] Richtmyer, R.D. and Morton, K.W. (1967), Difference Equations for Initial-Value Problems,
Interscience Publishers, New York.
[35] Young, D.L., “Time-dependent multi-material flow with large fluid distortion”, Numerical
Methods for Fluids Dynamics, Ed. K. W. Morton and M.J. Baines, Academic Press, New-York
(1982).
[36] Alia, A. and Souli M., 2006, “High explosive simulation using multi-material formulations”,
Applied Thermal Engineering, Vol. 26, pp.1032–1042.
[37] Woodward, P.R., and Collela, P., “The numerical simulation of two-dimensional fluid flow with
strong shocks”, Lawrence Livermore National Laboratory, UCRL-86952, (1982).
[38] Van Leer, B., “Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to
Numerical Convection”, Journal of Computational Physics 23 (1977) p 276–299.
[39] S. Rabier, M. Medale, Computation of free surface flows with a projection FEM in a moving
mesh framework, Comput. Methods Appl. Mech. Engrg. 192 (41–42) (2003) 4703–4721.
[40] D.J. Benson, An efficient, accurate, simple ALE method for nonlinear finite element programs,
Comput. Methods Appl. Mech. Engrg. 72 (1989) 305–350.
[41] Z.H. Zhong, Finite Element Procedures for Contact-impact Problems, Oxford University Press,
Oxford, 1993.
[42] T. Belytschko, M.O. Neal, Contact-impact by the pinball algorithm with penalty, projection, and
Lagrangian methods, in: Proceedings of the Symposium on Computational Techniques for
Impact, Penetration, and Performation of Solids AMD, vol. 103, ASME, New York, NY, 1989,
pp. 97–140.