Ale

Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

www.elsevier.com/locate/cma

ALE formulation for ¯uid±structure interaction problems M. Souli a,1, A. Ouahsine b,*, L. Lewin c b

a Livermore Software Technology Corp., Livermore, CA, USA 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 independently from material geometry. In ¯uid±structure interaction problems, where the ¯uid mesh near the structure undergoes large deformations 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 geometrical calculations associated with handling an arbitrary mesh. In this paper we apply the ALE concept to ¯uid±structure 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 calculations to continue. A most recent and interesting application of the ALE methods occurs in the medical ®eld, where experimental studies are very dicult to carry out, performing the experiments and extracting quantitative *

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

660

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

¯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 elastic±plastic 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 soil±structure interaction, detailed description of the method is given in [4]. Di€erent 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 algorithms, numerical results are in good agreement with experimental results, see [2]. For the 3D applications, 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 distorted 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 di€erential equations governing the mesh motion. In Section 3 di€erent 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 e€ectiveness 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 developed 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 dicult 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 con®guration time derivative, of …Xi ; t† of …xi ; t† of …xi ; t† ˆ ÿ wi ; ot ot oxi

…1†

where Xi is the Lagrangian coordinate, i the referential coordinate, xi the Eulerian coordinate, ui and wi are the material and reference velocities, respectively. Let denote by v the velocity of the material and by u the

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

661

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 ov oq ÿ wi : ˆ ÿq ot oxi oxi

…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 Navier±Stokes equations ov ovi ˆ ÿ…rij;j ‡ qbi † ÿ qwi : oxj ot

…3†

The stress tensor rij is described as follows rij ˆ ÿpdij ‡ l…vi;j ‡ vj;i †:

…4†

The last equations are solved with the following boundary conditions and initial conditions vi ˆ U 0

on C1 ;

rij ni ˆ 0

on C2 ;

…5†

where C1 [ C2 ˆ C;

C1 \ C2 ˆ 0;

…6†

C is the whole boundary of the calculation domain, and C1 and C2 are partial boundaries of C. The superscript means prescribed value, ni is the outward unit normal vector on the boundary, and dij is Kronecker's delta function. The velocity ®eld is assumed as known at t ˆ 0 in the whole domain X. vi …x; 0† ˆ 0:

…7†

(iii) The energy equation oE oE ˆ ÿ…rij vi;j ‡ qbi vj † ÿ qwj : ot oxj

…8†

Note that the Eulerian equations are derived by assuming that the velocity of the reference con®guration is zero and that the relative velocity between the material and the reference con®guration 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 dicult 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 di€erent 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: ovi ˆ rij;j ‡ qbi ; ot oE ˆ rij vi;j ‡ qbi vi : ot

…9† …10†

662

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

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 dicult form a standpoint of implementation, and is problem dependant. Di€erent 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 de®ne curvilinear coordinates n ˆ …n1 ; n2 ; n3 †, is given by r2 n ˆ 0:

…11†

We solve Eq. (11) for the coordinates x…ni †; …i ˆ 1; 2; 3† of the mesh lines: that is we invert them so that the geometric coordinates x ˆ …x1 ; x2 ; x3 † become the dependent variables and the curvilinear coordinates x…ni † …i ˆ 1; 2; 3† the independent variables. By the usual methods of changing variables we obtain: a1 on1 n1 x ‡ a2 on2 n2 x ‡ a3 on3 n3 x ‡ 2b1 on1 n2 x ‡ 2b2 on1 n3 x ‡ 2b3 on2 n3 x ˆ 0;

…12†

where ai ˆ oni x21 ‡ oni x22 ‡ oni x23 ;

i ˆ 1; 2; 3;

…13†

b1 ˆ …on1 x
on3 x†…on2 x
on3 x† ÿ …on1 x
on2 x†on3 x2 ;

…14†

b2 ˆ …on2 x
on1 x†…on3 x
on1 x† ÿ …on2 x
on3 x†on1 x2 ;

…15†

b3 ˆ …on3 x
on2 x†…on1 x
on2 x† ÿ …on3 x
on1 x†on2 x2 :

…16†

The equipotential method can be regarded as a minimization of the mesh density gradient by using variational 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 equipotential 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 successfully 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

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

n‡1 xSA ˆ

N 1 X xn : N iˆ1 i

663

…17†

3.3. Kikuchi's algorithm The volume weighting algorithm proposed by Kikushi, uses a volume weighted average of the coordinates of the centroids of the elements surrounding a node. First, we de®ne the coordinates of the element center xna ˆ

N 1X xn : 8 iˆ1 i

The new location of the node is given by PM Va xna n‡1 : xK ˆ Paˆ1 M aˆ1 Va

…18†

…19†

3.4. Combining smoothing algorithms The previous three algorithms de®ned above can be used simultaniously with di€erent scale factors. Let's denote by wE ; wSA and wK 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: n‡1 ‡ wSA xn‡1 xn‡1 ˆ wE xn‡1 E SA ‡ wK xK :

…20†

4. Advection process and Euler formulation One of the key points in ¯uid±structure 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 common 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 accompanied 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 di€usive than the simpler schemes [6,10,11]. The numerical approximation of the advective process requires important choices and compromises to minimize both arti®cial numerical di€usion and dispersion [13,18,22]. To illustrate these methods we consider, for simplicity, the one-dimensional scalar advection equation: o/ ou/ ‡ ˆ 0: ot ox

…21†

The discretization of Eq. (21) is obtained in the ¯ux form on a staggered grid by integration over the ®nite volume Xi ˆ ‰xiÿ12 ; xi‡12 Š and can be expressed explicitly in time as ˆ /ni ÿ /n‡1 i

i i Dt h Dt h n n …u/†ni‡1 ÿ …u/†niÿ1 ˆ /ni ÿ Fi‡1 ÿ Fiÿ 1 ; 2 2 2 2 Dx Dx

…22†

664

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

where F is the ¯ux of the scalar /, Dx the regular grid spacing , Dt the time step, n the time level tn ˆ nDt, n n /ni ˆ /…xi ; tn † and Fi‡ 1 , Fiÿ1 are the / numerical ¯uxes through the right and left boundaries of the grid cell, 2 2 respectively. It is known that the second order schemes eliminates the numerical di€usion, but they give rise to a nonphysical 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 di€erences, one can use the most straightforward approximation: /ni‡1 ˆ 0:5…/ni ‡ /ni‡1 †: 2

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: ( n n   ui‡1 /i if uni‡1 P 0; n 2 2 ˆ …23† Fi‡1 uni‡1 /ni‡1 if uni‡1 < 0; 2 Up 2

2

which can be written by a single formula      i 1 h n n /i‡1 ‡ uni‡1 /ni ‡ uni‡1 ÿ uni‡1 /ni‡1 : ˆ Fi‡ 1 2 2 Up 2 2 2 2

…24†

Since this approximation avoids non-physical oscillations, it develops the excessive numerical di€usion. The best way to overcome non-physical oscillations and excessive numerical di€usion, 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 /ni‡1 as the sum of the di€usive ®rst order upwind term and an ``anti-di€usive'' one. The higher order 2 anti-di€usive part is multiplied by the ¯ux-limiter, which depends locally on the nature of the solution by means of the non-linear function ri‡12 . The ¯ux-limiter scheme with the hybrid procedure is expressed as: 8 ÿ
 ‡  > < /ni ‡ 12 /ni‡1 ÿ /ni W ri‡ if ui‡12 P 0; 1 2   …25† /ni‡1 ˆ 2 > /n ÿ 1 ÿ/n ÿ /n
W rÿ : if u 1 < 0: 1 i‡ i‡1 i‡1 i 2 i‡ 2 2

The function ri‡12 refers to the slopes ratio at the neighborhood of the interfaces in the upwind direction (23). 8 n /i ÿ /niÿ1 > ‡ > > ˆ ri‡ if ui‡12 P 0; 1 > < /ni‡1 ÿ /ni 2 …26† ri‡12 ˆ n n > / ÿ / > i‡2 i‡1 ÿ > > ˆ ri‡1 if ui‡12 < 0: : n 2 /i‡1 ÿ /ni The interface value /niÿ1 may be deduced from /ni‡1 , by replacing the indices i by i ÿ 1. Note that from (25), 2 2 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 n‡1 Š 6 TV‰C n Š; X j di‡12 C n j; TV‰C n Š ˆ i

d1‡12 …:† ˆ …:†i‡1 ÿ …:†i :

…27†

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

665

Fig. 1. Numerical results of di€erent schemes compared with exact solution.

This leads to the restricting range: 06

W…r† 62 r

and

0 6 W…r† 6 2:

…28†

The various limiters schemes used in a brief numerical comparison, shown in Fig. 1, are: Van Leer [21]: /…r† ˆ …r‡ j r j†=…1‡ j r j†. 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 elastic±plastic 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 e€ects 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

666

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

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.

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 algorithm 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

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

667

Fig. 4. ALE formulation for Taylor bar impact, t ˆ 1000 ls:

Fig. 5. Comparison between Lagrangian and ALE time step.

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 ¯uid±structure 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 de®nes the pressure as a function of density and internal energy per unit volume, input parameters of this equation are given by Dobratz [8], for a variety of

668

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

Fig. 6. Fluid and structure meshes. Spherical pressure wave approaching the plate, t ˆ 50 ls:

Fig. 7. Fluid±structure interaction. Spherical pressure wave reaching the plate, t ˆ 200 ls:

high explosive materials. Away from the structure, an Eulerian multi-material formulation is used for the HE and the water, this formulation allows two di€erent 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 formulation. 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.

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

669

Fig. 8. Fluid±structure interaction. Plate deformation under pressure loading, with re¯ected pressure.

Fig. 9. Fluid±structure interaction. Initial and ®nal shapes of the plate.

5.2. Industrial example This example is directed toward the development of techniques for simulating fuel slosh in nonrectangular 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 optimizing the design of fuel slosh mitigating devices such as ba‚es. The central objective of these numerical simulations is to quantify the pitching moment due to the combined e€ects of pitch inertia and fuel sloshing. Fig. 10 shows the overall analysis model. The tank is

670

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

Fig. 10. ALE ¯uid mesh at t ˆ 500 ls:

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 approximately 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. E€ects 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 equilibrium. 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 becomes zero at t ˆ 3:4 s, after which the di€erence between the moment in the torsion spring and the static equilibrium value is due to the combined e€ects of pitch inertia and fuel slosh, neglecting the small moment due to pitch damping. An estimate of the incremental moment due the fuel slosh Ms , can be obtained by considering the pitch equation of motion, written as

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

671

Fig. 11. Analysis approach and ¯uid model. The tank is approximately half full of fuel.

 ‡ CH _ ‡ KH ˆ W d ‡ Ms ; IH

…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). Ms ˆ 57 483 ÿ 31 490 ÿ …398†…87:8 ÿ 30:03† ˆ 3000 in: lb:

…30†

As a check, an upper bound for Ms 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, Ms < …398†…30† ˆ 11 850 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 de®ne the

672

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

Fig. 12. Vertical force response during settling phase, ptop ˆ 1 psi:

Fig. 13. Pitch moment response during settling phase, ptop ˆ 1 psi:

¯uid±structure 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 ba‚es, and immersed object of varying shapes.

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

673

Fig. 14. Applied moment to initiate fuel slosh.

Fig. 15. Pitch moment history.

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.

674

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

Fig. 16. Inertial moment history.

3. The techniques developed can be used to optimize the use of slosh-mitigating devises such as ba‚es 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. Di€erent 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. Di€erent 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 ¯uid±structure 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 ecient and for the advection to be more accurate. In this paper, we show that the ALE methods we presented can be applied for ¯uid±structure 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 formulation 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

M. Souli et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 659±675

675

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) 235±294. [2] D.J. Benson, An ecient, accurate ALE method for nonlinear ®nite element programs, J. Comput. Methods Appl. Mech. Engrg. 72 (1989) 305±350. [3] D.J. Benson, A new two-dimensional ¯ux-limited shock viscosity, J. Comput. Methods Appl. Mech. Engrg. 93 (1991) 39±95. [4] I. Shahrour, B. Benchekh, Analysis of the soil±structure 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) 329±349. [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) 217±284. [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) 155±196. [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) 155±196. [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) 285±312. [11] T. Nomura, T.J.R. Hughes, An arbitrary Lagrangian±Eulerian ®nite element method for interaction of ¯uid and a rigid body, J. Comput. Methods Appl. Mech. Engrg. 95 (1992) 115±138. [12] V. Palanisamy, M. Kawahara, A fractional step arbitrary Lagrangian Eulerian ®nite element, Comp. Fluid Dynamics 1 (1993) 57±77. [13] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Engng. 19 (1979) 59±98. [14] P.K. Sweby, High-resolution schemes using ¯ux-limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (1984) 995±1011. [15] A. Harten, High-resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357±393. [16] M. Souli, J.P. Zolesio, Shape derivative of descretized problems, Comput. Methods Appl. Mech. Engrg. 108 (1993) 187±199. [17] M. Souli, J.P. Zolesio, Finite element method for free surface ¯ow problems, Comput. Methods Appl. Mech. Engrg. 129 (1996) 43±51. [18] R. Leveque, High-resolution conservative algorithms for advection in incompressible ¯ow, SIAM J. Numer. Anal. 33 (1996) 627±665. [19] B. Roe, Some contributions to the modeling of discontinuous ¯ows, Lectures in Appl. Math. 22 (1985) 163±193. [20] B. Roe, D. Sidilkover, Optimum positive linear schemes for advection in two and three dimensions, SIAM J. Numer. Anal. 29 (1992) 1542±1568. [21] B. Van Leer, Towards the ultimate conservative di€erence scheme. II Monotonicity and conservation combined in a second order scheme, J. Comput. Phys. 14 (1974) 361±370. [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) 76±84.