ALE formulation for uidstructure interaction problems

M. Souli
a,1
, A. Ouahsine
b,
*
, L. Lewin
c
a
Livermore Software Technology Corp., Livermore, CA, USA
b
Laboratoire de M ecanique de Lille, Univ. des Scienc. & Techn. de Lille, URA CNRS 1441, Bd. Paul Langevin,
59655 Villeneuve d'Ascq Cedex, France
c
The Boeing Co., Aircraft Division, Seattle, Washington, USA
Received 21 April 1999
Abstract
Arbitrary Lagrangian Eulerian (ALE) nite element methods gain interest for the capability to control mesh geometry indepen-
dently from material geometry. In uidstructure interaction problems, where the uid mesh near the structure undergoes large de-
formations and becomes unacceptably distorded, which drive the time step to a very small value for explicit calculations, the ALE
methods or rezoning are used to create a new undistorted mesh for the uid domain, which allows the calculations to continue. The
mathematical basis of the ALE and rezoning algorithms is simple, but their implementation is complicated due to the tedious geo-
metrical calculations associated with handling an arbitrary mesh. In this paper we apply the ALE concept to uidstructure interaction
problems. We will explain the underlying ideas of the method and a possible way to control the distortion of the mesh is given. Results
of an academic as well as an industrial problem are presented. 2000 Elsevier Science S.A. All rights reserved.
1. Introduction
The arbitrary Lagrangian Eulerian (ALE) approach is based on the arbitrary movement of a reference
domain which, additionally to the common material domain and spatial domain, is introduced as a third
domain, as detailed in [5]. In this reference domain, which will later on correspond to the nite 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 and large deformations, and interface contact problems.
The purpose of the paper is to describe new ALE algorithms that have been developed and used in
hydrocodes; the original motivation for developing these codes was for solving defense problems, a detailed
description of the motivation is given by Benson [1]. The range of applications has been exented in recent
years. Current applications include sloshing tank problems involving free surfaces and high velocity impact
problems, where the target is treated as a uid material which undergoes large deformations, and the
penetrator is structure material. For these problems, a contact with eroding elements is used at the interface
separating the target from the penetrator, and ALE algorithms are used on the uid mesh near the structure
to prevent the mesh to get severely distorted and the time step from decreasing which prevent the calcu-
lations to continue.
A most recent and interesting application of the ALE methods occurs in the medical eld, where ex-
perimental studies are very dicult to carry out, performing the experiments and extracting quantitative
www.elsevier.com/locate/cma
Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
*
Corresponding author Tel.: +33-03-20337152; fax: +33-03-20337153.
E-mail address: [email protected] (A. Ouahsine).
1
Present address: Universite d'Artois, LAMH-FSA Bethune, 62400 Bethune, France.
0045-7825/00/$ - see front matter 2000 Elsevier Science S.A. All rights reserved.
PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 4 3 2 - 6
ow data is time consuming and expensive process. Examples are the blood vessels pulsingly own by the
blood [9]. In these applications the blood is modeled as a uid using a non-Newtonian material model, and
the vessels as a structure using an elasticplastic material model. For high pressure, the vessels undergo
large deformations which leads to the remesh of the uid domain.
Other applications of this formulation were used for soilstructure interaction, detailed description of
the method is given in [4]. Dierent types of contact have been used to simulate the interface between the
uid and the structure; the uid mesh close to the contact interface is treated using an ALE formulation, in
order to adapt the uid mesh to the structure mesh for the contact algorithms. while away from the contact
an Eulerian formulation on a xed mesh for the uid is suitable.
For 2D applications where the remesh of the uid domain can be easily treated with actual ALE al-
gorithms, numerical results are in good agreement with experimental results, see [2]. For the 3D applica-
tions, the remesh of the uid domain is more complex due to the complexity of the 3D geometry of
computational domain. The actual ALE algorithms are not general and robust enough to handle complex
3D meshes, and fail for most complex 3D problems. For most of the 3D applications, automatic mesh
generators are called internally to create a new mesh with a new topology, given the boundary segments,
this method refers to a rezoning method, the dependent variables: velocity, pressure, internal energy, stress
components and plastic strain are updated on the new mesh by using a remap algorithm. Unlike a rezoning,
the topology of the mesh is xed in an ALE method. The accuracy of an ALE calculation is often superior
to the accuracy of a rezoned calculation because the algorithms used to remap the solution from the dis-
torted to the undistorted mesh is second order accurate for the ALE formulation when using second order
advection algorithms, while the algorithm for the remap in the rezoning is only rst order accurate.
The outline of this paper is arranged as follows. In Section 2, a general description of the ALE method is
introduced, we describe the partial dierential equations governing the mesh motion. In Section 3 dierent
ALE algorthims are described, some of them are classical and can be found in published papers [1,2,12], the
other ones not published, and based on geometrical basis, developed at Livermore Software Technology
Corporation, in order to solve practical problems. An application of the algorithm is described in detail
that indicates the eectiveness of the approach. Section 4 is devoted to rst and second advection schemes.
Numerical results are compared to analytical results. In Section 5, two academic numerical examples and
one industrial application are presented. The rst academic examples concern, the classical Taylor problem
(the bar impact), and the second concern the uid structure interaction problem involving underwater
explosion. In the industrial application, the simulation of fuel slosh in aircraft thanks gas has been de-
veloped and illustrated for a tank with curved surfaces and approximately half full of fuel.
2. Governing processes
An ALE formulation contains both pure Lagrangian and pure Eulerian formulations. The pure
Lagrangian description is the approach that: the mesh moves with the material, making it easy to track
surfaces and apply boundary conditions. Using an Eulerian description, the mesh remains xed while the
material passes through it. Surfaces and boundary conditions are dicult to track using this approach;
however, mesh distortion is not a problem because the mesh never changes. In solid mechanics a pure
Eulerian formulation it is not useful because it can handle only a single material in an element, while an
ALE formulation is assumed to be capable of handling a single material in an element.
In the ALE description, an arbitrary referential coordinate is introduced in addition to the Lagrangian
and Eulerian coordinates. The material derivative with respect to the reference coordinate can be described
as Thus the ALE equations are derived by substituting the relationship between the material time derivative
and the reference conguration time derivative,
of (X
i
; t)
ot
=
of (x
i
; t)
ot
w
i
of (x
i
; t)
ox
i
; (1)
where X
i
is the Lagrangian coordinate, i the referential coordinate, x
i
the Eulerian coordinate, u
i
and w
i
are
the material and reference velocities, respectively. Let denote by v the velocity of the material and by u the
660 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
velocity of the mesh. In order to simplify the equations we introduce the relative velocity w = v u. Thus
the governing equations for the ALE formulation are given by:
(i) The conservation of mass equation.
oq
ot
= q
ov
ox
i
w
i
oq
ox
i
: (2)
(ii) The momentum equation. The strong form of the problem governing Newtonian uid ow in a xed
domain consists of the governing equations and suitable initial and boundary conditions. The equations
governing the uid problem are the ALE description of the NavierStokes equations
ov
ot
= (r
ij;j
qb
i
) qw
i
ov
i
ox
j
: (3)
The stress tensor r
ij
is described as follows
r
ij
= pd
ij
l(v
i;j
v
j;i
): (4)
The last equations are solved with the following boundary conditions and initial conditions
v
i
= U
0
on C
1
;
r
ij
n
i
= 0 on C
2
;
(5)
where
C
1
C
2
= C; C
1
C
2
= 0; (6)
C is the whole boundary of the calculation domain, and C
1
and C
2
are partial boundaries of C. The
superscript means prescribed value, n
i
is the outward unit normal vector on the boundary, and d
ij
is
Kronecker's delta function. The velocity eld is assumed as known at t = 0 in the whole domain X.
v
i
(x; 0) = 0: (7)
(iii) The energy equation
oE
ot
= (r
ij
v
i;j
qb
i
v
j
) qw
j
oE
ox
j
: (8)
Note that the Eulerian equations are derived by assuming that the velocity of the reference conguration
is zero and that the relative velocity between the material and the reference conguration is therefore the
material velocity. The term in the relative velocity is usually referred to as the advective term, and accounts
for the transport of material past the mesh. It is the additional term in the equations that makes solving the
ALE equations much more dicult numerically than the Lagrangian equations, where the relative velocity
is zero.
There are two ways to implement the ALE equations, and they correspond to the two approaches
taken in implementing the Eulerian viewpoint in uid mechanics. The rst way solves the fully
coupled equations for computational uid mechanics, this approach used by dierent authors can
handle only a single material in an element. The alternative approach is referred to as an operator
split in the literature, where the calculation, for each time step is divided into two phases. First a
Lagrangian phase is performed, in which the mesh moves with the material, in this phase the changes
in velocity and internal energy due to the internal and external forces are calculated. The equilibrium
equations are:
ov
i
ot
= r
ij;j
qb
i
; (9)
oE
ot
= r
ij
v
i;j
qb
i
v
i
: (10)
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 661
In the second phase, the advection phase, transport of mass, internal energy and momentum across cell
boundaries are computed; this may be thought of as remapping the displaced mesh at the Lagrangian phase
back to its original or arbitrary position.
3. Smoothing algorithms
Deciding where to move the nodes is a very challenging step, this step is the most dicult form a
standpoint of implementation, and is problem dependant. Dierent relaxation stencils have been used in
hydrocodes, the original one is the equipotential stencil proposed by Winslow [7] and implemented in
LSDYNA was found to be stable for a broad range of problems. Other ALE methods which are purely
geometrical, the simple average algorithm, and the weighted volume, the Kikucshi algorithm have been
intensively used in LSDYNA to solve underwater explosion problems and explosion air, these algorithms
are combined with the equipotential algorithm using scale factors for each algorithm.
3.1. Equipotential algorithm
In this algorithm, the inverse of a Laplace equation is solved. A stencil algorithm for solving the
Laplacian equation is derived by Winslow [7]. The boundary nodes move tangentially to the boundary
using a 2D equipotential smoothing algorithm. The boundary nodes are moved before the interior mesh
is smoothen. The inverted form of Laplacian equation, in three dimensions where we dene curvilinear
coordinates n = (n
1
; n
2
; n
3
), is given by
\
2
n = 0: (11)
We solve Eq. (11) for the coordinates x(n
i
); (i = 1; 2; 3) of the mesh lines: that is we invert them so that
the geometric coordinates x = (x
1
; x
2
; x
3
) become the dependent variables and the curvilinear coordinates
x(n
i
) (i = 1; 2; 3) the independent variables. By the usual methods of changing variables we obtain:
a
1
o
n
1
n
1
x a
2
o
n
2
n
2
x a
3
o
n
3
n
3
x 2b
1
o
n
1
n
2
x 2b
2
o
n
1
n
3
x 2b
3
o
n
2
n
3
x = 0; (12)
where
a
i
= o
n
i
x
2
1
o
n
i
x
2
2
o
n
i
x
2
3
; i = 1; 2; 3; (13)
b
1
= (o
n
1
x o
n
3
x)(o
n
2
x o
n
3
x) (o
n
1
x o
n
2
x)o
n
3
x
2
; (14)
b
2
= (o
n
2
x o
n
1
x)(o
n
3
x o
n
1
x) (o
n
2
x o
n
3
x)o
n
1
x
2
; (15)
b
3
= (o
n
3
x o
n
2
x)(o
n
1
x o
n
2
x) (o
n
3
x o
n
1
x)o
n
2
x
2
: (16)
The equipotential method can be regarded as a minimization of the mesh density gradient by using vari-
ational formulation methods to minimize a functional cost, see [16].
3.2. Simple average
A simple averaging of the location of the surrounding nodes allows the mesh to be more uniform, since
this algorithm is nonlinear, few iterations are processed to smooth the mesh. Combined with the equipo-
tential algorithm, with a smaller scale factor for the simple average, this algorithm allows the boundary
nodes to adapt to the new mesh created by the equipotential algorithm. This combination has been suc-
cessfully used in LSDYNA for a broad range of problems, including impact problems, and underwater
explosion problems using a single material ALE formulation. One of the inconveniences of using the simple
average, is that it tends to eliminate the mesh grading in the material domain, and tends to smooth out any
steep gradient between elements. The new location of the node is given by
662 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
x
n1
SA
=
1
N
X
N
i=1
x
n
i
: (17)
3.3. Kikuchi's algorithm
The volume weighting algorithm proposed by Kikushi, uses a volume weighted average of the coordi-
nates of the centroids of the elements surrounding a node. First, we dene the coordinates of the element
center
x
n
a
=
1
8
X
N
i=1
x
n
i
: (18)
The new location of the node is given by
x
n1
K
=
P
M
a=1
V
a
x
n
a
P
M
a=1
V
a
: (19)
3.4. Combining smoothing algorithms
The previous three algorithms dened above can be used simultaniously with dierent scale factors. Let's
denote by w
E
; w
SA
and w
K
the scale factors, respectively used for the equipotential, the simple average and
Kikuchi's algorithm. These factors are positive and smaller than one. The nale location of the node is
given by:
x
n1
= w
E
x
n1
E
w
SA
x
n1
SA
w
K
x
n1
K
: (20)
4. Advection process and Euler formulation
One of the key points in uidstructure interaction is the necessity to apply the ALE formulation. This
approach is based on the arbitrary displacement of the reference domain which, additionally to the com-
mon material domain and spatial domain, is introduced as a third domain. To the reference domain
corresponds the nite element mesh. Hence, the arbitrary movement of the reference frame must be ac-
companied by reliable mesh moving algorithms, capable to deal with convenient moving boundaries, free
surfaces [17] and large deformations. Thus when applying the ALE formulations, both geometrical as well
as advective processes have to be overcome.
In this section we investigate several ux-limiter schemes, that have been chosen because they satisfy
many of the requirements of a good advection scheme, they are total variation diminishing (TVD), mass
conservative and less diusive than the simpler schemes [6,10,11].
The numerical approximation of the advective process requires important choices and compromises to
minimize both articial numerical diusion and dispersion [13,18,22]. To illustrate these methods we
consider, for simplicity, the one-dimensional scalar advection equation:
o/
ot

ou/
ox
= 0: (21)
The discretization of Eq. (21) is obtained in the ux form on a staggered grid by integration over the
nite volume X
i
= [x
i
1
2
; x
i
1
2
[ and can be expressed explicitly in time as
/
n1
i
= /
n
i

Dt
Dx
u/ ( )
n
i
1
2
h
u/ ( )
n
i
1
2
i
= /
n
i

Dt
Dx
F
n
i
1
2
h
F
n
i
1
2
i
; (22)
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 663
where F is the ux of the scalar /, Dx the regular grid spacing , Dt the time step, n the time level t
n
= nDt,
/
n
i
= /(x
i
; t
n
) and F
n
i
1
2
, F
n
i
1
2
are the / numerical uxes through the right and left boundaries of the grid cell,
respectively.
It is known that the second order schemes eliminates the numerical diusion, but they give rise to a non-
physical oscillations near regions of large gradients. In the next subsection we shall present a method to
obtain higher order schemes without oscillations.
4.1. Flux-limiter methods
To approximate the partial derivative ou/=ox by a central dierences, one can use the most straight-
forward approximation:
/
n
i
1
2
= 0:5(/
n
i
/
n
i1
):
Unfortunately, this approximation gives rise to spurious oscillations that perturb the numerical solution.
These free-oscillations will be provided by using the following upstream scheme:
F
n
i
1
2

Up
=
u
n
i
1
2
/
n
i
if u
n
i
1
2
P0;
u
n
i
1
2
/
n
i1
if u
n
i
1
2
< 0;
(
(23)
which can be written by a single formula
F
n
i
1
2

Up
=
1
2
/
n
i
1
2
h
u
n
i
1
2

/
n
i
u
n
i
1
2

u
n
i
1
2

/
n
i1
i
: (24)
Since this approximation avoids non-physical oscillations, it develops the excessive numerical diusion.
The best way to overcome non-physical oscillations and excessive numerical diusion, is the use of an
hybrid scheme. It consists on using the second order numerical ux in smooth regions and the rst order
one (the monotonic upwind method) in the vicinity of discontinuities. This procedure is carried out by
introducing a ux-limiter schemes, based on the local gradient of the solution (23). We write the interface
value /
n
i
1
2
as the sum of the diusive rst order upwind term and an ``anti-diusive'' one. The higher order
anti-diusive part is multiplied by the ux-limiter, which depends locally on the nature of the solution by
means of the non-linear function r
i
1
2
. The ux-limiter scheme with the hybrid procedure is expressed as:
/
n
i
1
2
=
/
n
i

1
2
/
n
i1
/
n
i

W r

i
1
2

if u
i
1
2
P0;
/
n
i1

1
2
/
n
i1
/
n
i

W r

i
1
2

if u
i
1
2
< 0:
8
>
<
>
:
(25)
The function r
i
1
2
refers to the slopes ratio at the neighborhood of the interfaces in the upwind direction (23).
r
i
1
2
=
/
n
i
/
n
i1
/
n
i1
/
n
i
= r

i
1
2
if u
i
1
2
P0;
/
n
i2
/
n
i1
/
n
i1
/
n
i
= r

i
1
2
if u
i
1
2
< 0:
8
>
>
>
>
<
>
>
>
>
:
(26)
The interface value /
n
i
1
2
may be deduced from /
n
i
1
2
, by replacing the indices i by i 1. Note that from (25),
if W = 0 one can obtain the upwind scheme, while if W = 1 the scheme is reduced to the centered one. For
stability, these ux-limiter schemes must satisfy the TVD concept (total variation diminishing), outlined in
[14,15,18],
TV[C
n1
[ 6TV[C
n
[;
TV[C
n
[ =
X
i
[ d
i
1
2
C
n
[;
d
1
1
2
(:) = (:)
i1
(:)
i
:
(27)
664 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
This leads to the restricting range:
0 6
W(r)
r
62 and 0 6W(r) 62: (28)
The various limiters schemes used in a brief numerical comparison, shown in Fig. 1, are:
Van Leer [21]: /(r) = (r [ r [)=(1 [ r [).
Superbee [19,20]: /(r) = max(0; min(1; 2r); min(2; r)).
5. Numerical examples and results
5.1. Academic examples
5.1.1. Example 1: The Taylor bar impact problem
The Taylor bar impact problem is a benchmark problem that appears frequently in the literature [3]. The
problem consists of a cylinder impacting a rigid wall (Fig. 2) with the impact velocity is 10 m/s. Since the
problem is axisymmetric, only one-fourth of the cylinder is modeled; a sliding contact interface is imposed
between the cylinder and the rigid wall. Eighteen hundred elements are used to model one-fourth of the
cylinder with a radius of 40 mm and a length of 100 mm. An elasticplastic hydrodynamic material with a
Gruneisen equation of state is used; with a shear modulus of 0.1725 Mb and a yield stress of 0:3410
3
Mb.
The problem is used to examine the eects of the ALE smoothing algorithms on the mesh, as well as the
time step. In explicit calculations, the time step is based on the characteristic of the elements, for the
Lagrangian calculations, the elements close to contact interface undergo large deformation, these elements
stretch in length and shrink in width and the size of the time step decreases, when the time step becomes too
Fig. 1. Numerical results of dierent schemes compared with exact solution.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 665
small, continuing the analysis becomes prohibitely expensive. A combination of the equipotential and
simple average algorithms with scale factors of 0.9 and 0.1 are used for this example. The smoothing al-
gorithm is performed every 10 cycles; these values are a reasonable set and are not to be taken as optimal.
Fig. 2 shows the initial mesh of one-fourth of the cylinder and Figs. 3 and 4 show the deformed mesh of
the Lagrangian and ALE calculations, after 100 ls. For this example, 8736 explicit cycles are needed in the
Lagrangian calculation (Fig. 3) to achieve 1000 ls for termination time. The time step starts with a value of
1.28 ls and drops to a value of 7:4810
2
ls. While the ALE calculation (Fig. 4) requires only 3648 cycles to
achieve 1000 ls. The time step (Fig. 5) starts with a value of 1.28 ls, and drops to 0.22 ls, which is about
three times higher than the Lagrangian time step (Fig. 3). The Lagrangian mesh undergoes more distortions
Fig. 2. Taylor bar impact at initial time. Model of one-fourth of the cylinder.
Fig. 3. Lagrangian formulation for Taylor bar impact, t = 1000 ls.
666 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
than the ALE mesh, and takes 22 mN on DEC Alpha 3000 workstation, whereas the ALE calculations
takes only 11 mN.
5.1.2. Example 2: The uidstructure coupling
This example (Fig. 6) illustrates an underwater explosion, using a C4 hgh eplosive (HE) to generate a
pressure wave after detonation. A chemical energy stored in the HE, is transformed into mechanical energy,
which generates a spherical pressure wave moving through the water medium and loading a thick structure
plate (Figs. 7 and 8). The plate is represented by several elements through the thickness to adequately
represent the bending. In this example the HE is modeled using burn program with a Jones, Wilkins, Lee
(JWL) equation of state. The JWL equation of state denes the pressure as a function of density and in-
ternal energy per unit volume, input parameters of this equation are given by Dobratz [8], for a variety of
Fig. 4. ALE formulation for Taylor bar impact, t = 1000 ls:
Fig. 5. Comparison between Lagrangian and ALE time step.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 667
high explosive materials. Away from the structure, an Eulerian multi-material formulation is used for the
HE and the water, this formulation allows two dierent materials within a single element. Close to the
structure the uid nodes at the interface move along the structure and have the same normal velocity as
the structure nodes (Figs. 8 and 9). This interface condition is enforced by the contact between the uid and
the structure. The uid is treated using an ALE formulation, and the structure with a Lagrangian for-
mulation. Fig. 10 illustrates the mesh structure after 500 ls, the deformed ALE mesh still keeps a good
shape for the calculations to continue further in time if desired.
Fig. 6. Fluid and structure meshes. Spherical pressure wave approaching the plate, t = 50 ls:
Fig. 7. Fluidstructure interaction. Spherical pressure wave reaching the plate, t = 200 ls:
668 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
5.2. Industrial example
This example is directed toward the development of techniques for simulating fuel slosh in non-
rectangular tank having curved boundaries and approximately half full of uid. Fuel slosh can lead to
undesirable pitch oscillations and result in instability. The simulation of this phenomena can aid in opti-
mizing the design of fuel slosh mitigating devices such as baes.
The central objective of these numerical simulations is to quantify the pitching moment due to the
combined eects of pitch inertia and fuel sloshing. Fig. 10 shows the overall analysis model. The tank is
Fig. 8. Fluidstructure interaction. Plate deformation under pressure loading, with reected pressure.
Fig. 9. Fluidstructure interaction. Initial and nal shapes of the plate.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 669
assumed rigid, since structural deformation can be considered negligible, and is given the overall pitch
inertia properties of the vehicle whose center of mass is 42 in. aft of the fuel tank. In the present study the
pressure on the upper surface of the uid is limited to 1 psi to preserve numerical stability. An equal
pressure is applied to the bottom of the tank. A grounded torsion spring at the center of mass is included to
provide a vehicle pitch frequency of 0.75 Hz and also provides a convenient means of measuring pitching
moment. Vehicle motion is limited to pitch and vertical translation. The tank is assumed to be approxi-
mately half full of fuel, shown in Fig. 10. The uid constitutive model employs a Gruneisen equation of
state. Material parameters are taken as those of water. Fluid motion can be three-dimensional and is
limited only by the shape of the tank. Eects of external aerodynamic and attitude control system forces are
not included in the analysis. The simulation is conducted in two phases, as described below.
5.2.1. Phase 1: Settling
In this rst phase of the simulation both gravity and pressure loads are applied as quasi-static functions
of time. Pressure is also applied to the bottom of the base of the tank to avoid unrealistically large pitch
motions during the settling phase. At the completion of settling, the system is in a state of static equilib-
rium. The load in the vertical spring at the center of mass should be equal to the weighing of the uid, 391.5
lbs, and the moment in the torsion spring should be roughly equal to the product of this load and the
distance from the center of mass of the uid to the axis of rotation; this moment is computed to be 22 617
in. lb. As shown in Figs. 11 and 12, which show net vertical force and pitching moment during settling, the
net vertical force agrees closely with the uid weight and the moment is about 15% above the computed
value. This is attributed primarily to the fact that the bottom of the tank is not at and the resultant force
on this surface has a horizontal as well as vertical component.
5.2.2. Phase 2: Fuel slosh
It is assumed that fuel slosh is initiated by a pitching moment imposed by an external aerodynamic
disturbance and/or the vehicle attitude control system. Here the moment is applied as a half-sine impulse,
shown in Fig. 13; its magnitude is set to produce an equivalent static pitch rotation of three degrees. During
Phase 2 pitch damping is reduced from 10% to 1% of critical. Fig. 14 shows the time variation of pitching
moment about the vehicle center of gravity as measured by the torsion spring. The applied moment be-
comes zero at t = 3:4 s, after which the dierence between the moment in the torsion spring and the static
equilibrium value is due to the combined eects of pitch inertia and fuel slosh, neglecting the small moment
due to pitch damping. An estimate of the incremental moment due the fuel slosh M
s
, can be obtained by
considering the pitch equation of motion, written as
Fig. 10. ALE uid mesh at t = 500 ls:
670 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
I

H C
_
H KH = Wd M
s
; (29)
where W and d are the uid weight and distance from the uid center of gravity, to the center of rotation,
respectively, and I is the pitch inertia. The variation of I due to uid motion is very small. Using the data of
Figs. 14 and 15, the latter showing the inertial moment, gives for the moment due to fuel slosh at a time of
4.2 s (see Fig. 16).
M
s
= 57483 31490 (398)(87:8 30:03) = 3000 in: lb: (30)
As a check, an upper bound for M
s
can be established by assuming that all the uid moves to the end of
the tank furthest away from the axis of rotation, a distance of approximately 30 in. Hence,
M
s
< (398)(30) = 11850 in: lb: (31)
Thus the above value is about 26% of the upper bound.
5.2.3. Discussion of results
In this numerical simulation, the uid is modeled with brick elements and the tank with shell elements,
using rigid body material. Fluid nodes in contact with tank are common nodes. Since the uid motion
along the tank is relatively small, common nodes between the uid and the tank can be used to dene the
Fig. 11. Analysis approach and uid model. The tank is approximately half full of fuel.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 671
uidstructure interface. In this example, the equipotential algorithm is used every 10 cycles in the analysis,
to allow the uid mesh close to the structure to move smoothly, and prevent time step from decreasing.
During the entire analysis the tine step is constant, and is controlled by the uid elements, since the tank is a
rigid body material. This analysis has been used to optimize devices in sloshing tanks, such as baes, and
immersed object of varying shapes.
Fig. 13. Pitch moment response during settling phase, p
top
= 1 psi:
Fig. 12. Vertical force response during settling phase, p
top
= 1 psi:
672 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
The analysis described above leads to the following conclusions:
1. An analytical technique for simulation of fuel slosh in aircraft thanks gas been developed and illustrated
for a tank with curved surfaces and approximately half full of fuel.
2. The simulation can be used for thanks of a wide variety of shapes and is well-suited for investigation of
the pitching moment attributable to fuel slosh.
Fig. 15. Pitch moment history.
Fig. 14. Applied moment to initiate fuel slosh.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 673
3. The techniques developed can be used to optimize the use of slosh-mitigating devises such as baes and
to simulate fuel slosh in tanks containing immersed objects of varying shapes.
6. Conclusion
The purpose of the ALE formulation presented in this paper is to demonstrate the ability of the ALE
smoothing method to continue the calculation beyond the point where the pure Lagrangian calculation
fails.
For most problems, The ALE formulation reduces the cost of an analysis involving large deformations.
Dierent smoothing approaches are adopted in this paper in order to remesh the computational domain
without requiring the user intervention. This also allows to reduce manual rezoning for problems, where
manual rezoning is necessary. Dierent smoothing algorithms are described. The equipotential methods
which is commonly used in the industry and academia to remap a distorted mesh into an undistorted mesh.
The simple average and the volume weighted algorithms are not very common and have been successfully
used in LSDYNA to solve problems involving uidstructure interaction. In order for these algorithms to
be used successfully, one should not allow the mesh to be highly distorted. Remap algorithms need to be
evoked more frequently for the smoothing algorithm to be more ecient and for the advection to be more
accurate.
In this paper, we show that the ALE methods we presented can be applied for uidstructure problems
coming out of academia and industrial application.
For the Taylor bar problem, the ALE calculations can continue beyond termination time if desired,
while the Lagrangian formulation shows that continuing the calculation will be prohibitely expensive, and
sometimes impossible.
The blast problem illustrated in the second example, shows the necessity of applying an Eulerian for-
mulation away the structure and ALE formulation closed to the structure. Since the mesh is rectangular
and the pressure waves is spherical, which mean that the nodes moves radially in rectangular mesh. The
Fig. 16. Inertial moment history.
674 M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675
Lagrangian formulation cannot handle the problem more than few cycles before the element get inverted
and the calculation stopped.
References
[1] D.J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, J. Comput. Methods Appl. Mech. Engrg. 99 (1992)
235294.
[2] D.J. Benson, An ecient, accurate ALE method for nonlinear nite element programs, J. Comput. Methods Appl. Mech. Engrg.
72 (1989) 305350.
[3] D.J. Benson, A new two-dimensional ux-limited shock viscosity, J. Comput. Methods Appl. Mech. Engrg. 93 (1991) 3995.
[4] I. Shahrour, B. Benchekh, Analysis of the soilstructure interaction under monotonic and cyclic loadings, in: Ch. Hirsh, O.C.
Zienkieviwicz, E. Onate (Eds.), Proceedings of the First European Conference on Numerical Methods in Engineering, Bruxelles,
Elsevier, Amsterdam, 1992.
[5] T.J.R. Hughes, W.K. Liu, T.K. Zimmerman, Lagrangian Eulerian nite elements formulation for viscous ows, J. Comput.
Methods Appl. Mech. Engrg. 21 (1981) 329349.
[6] T.J.R. Hughes, T.E. Tezduar, Finite element methods for rst-order hyperbolic systems with particular emphasis on the
compressible Euler equations, J. Comput. Methods Appl. Mech. Engrg. 45 (1984) 217284.
[7] A.M. Winslow, Equipotential zoning of the interior of a three-dimensional mesh, Lawrence Radiation Laboratory, UCRL-7312,
1990.
[8] B.M. Dobratz, Lawrence Livermore National Laboratory, Explosives Handbook. Properties of Chemical Explosives and
Explosive Stimulants, University of California, LLNL Report ICRL-52997, 1981; J. Comput. Methods Appl. Mech. Engrg. 158
(1981) 155196.
[9] Ch.A. Taylor, T.J.R. Hughes, C.K. Zarins, Finite element modeling of blood ow in arteries, J. Comput. Methods Appl. Mech.
Engrg. 158 (1998) 155196.
[10] T.J.R. Hughes, M. Mallet, M. Mizukami, A new nite element formulation for computational uid dynamics: II. Beyond SUPG,
J. Comput. Methods Appl. Mech. Engrg. 45 (1984) 285312.
[11] T. Nomura, T.J.R. Hughes, An arbitrary LagrangianEulerian nite element method for interaction of uid and a rigid body,
J. Comput. Methods Appl. Mech. Engrg. 95 (1992) 115138.
[12] V. Palanisamy, M. Kawahara, A fractional step arbitrary Lagrangian Eulerian nite element, Comp. Fluid Dynamics 1 (1993)
5777.
[13] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput.
Methods Appl. Mech. Engng. 19 (1979) 5998.
[14] P.K. Sweby, High-resolution schemes using ux-limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (1984)
9951011.
[15] A. Harten, High-resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357393.
[16] M. Souli, J.P. Zolesio, Shape derivative of descretized problems, Comput. Methods Appl. Mech. Engrg. 108 (1993) 187199.
[17] M. Souli, J.P. Zolesio, Finite element method for free surface ow problems, Comput. Methods Appl. Mech. Engrg. 129 (1996)
4351.
[18] R. Leveque, High-resolution conservative algorithms for advection in incompressible ow, SIAM J. Numer. Anal. 33 (1996)
627665.
[19] B. Roe, Some contributions to the modeling of discontinuous ows, Lectures in Appl. Math. 22 (1985) 163193.
[20] B. Roe, D. Sidilkover, Optimum positive linear schemes for advection in two and three dimensions, SIAM J. Numer. Anal. 29
(1992) 15421568.
[21] B. Van Leer, Towards the ultimate conservative dierence scheme. II Monotonicity and conservation combined in a second order
scheme, J. Comput. Phys. 14 (1974) 361370.
[22] G.D. Van Albada, B. Van Leer, J.W.W. Roberts, A comparative study of computational methods in cosmic gas dynamics,
Astronom. Astrophys. 108 (1982) 7684.
M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659675 675