- Email: [email protected]
Arbitrary Lagrangian–Eulerian and free surface methods in fluid mechanics
Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
www.elsevier.com/locate/cma
Arbitrary Lagrangian±Eulerian and free surface methods in ¯uid mechanics M. Souli a
a,*
, J.P. Zolesio
b
Lab. de Mechanique, Universit e des Sciences et Technologie de Lille 1, Bd Paul Langevin, Cite Scienti®que, 59655 Lille, France b Ecole des Mines de Paris, Sophia Antipolis, Route des Lucioles, Valbonne, France Received 25 March 2000; received in revised form 25 August 2000; accepted 30 May 2001
Abstract A computational procedure is developed to solve problems of viscous incompressible ¯ows under large free surface motions. The arbitrary Lagrangian±Eulerian (ALE) method is used to move the free surface nodes as well as the internal nodes. The coupling of the mesh motion equations and the ¯uid equations is essentially done through the free surface boundary conditions. An elastic model with free surface boundary conditions is used to model the mesh motion. The boundary condition at the free surface is deduced from the fact that no ¯uid particle can cross the free surface. The ¯uid problem is solved using a Galerkin/least-square (GLS) ®nite element formulation. Ó 2001 Elsevier Science B.V. All rights reserved.
1. Introduction In continuum mechanics, two descriptions are considered for the motion in a continuum media: The ®rst is the Eulerian description, where we focus attention on a particular volume in space. The volume is ®xed with respect to a laboratory frame, and we study the ¯uid as it passes through the ®xed volume. The description is one in which the ¯uid is continuously renewed inside the volume, the Eulerian description is not the simplest in which the basic equations of ¯uid motion can be formulated. A convective term is introduced to express the material time derivative in the reference con®guration. The convective term gives a nonsymmetrical form of the Galerkin formulation. Since the computational domain is ®xed, the Eulerian description has the advantage of preserving the mesh regularity. The second is the Lagrangian description, in which we identify and follow a particular region of ¯uid. The volume of ¯uid changes in shape, while the total mass remains constant. In the Lagrangian description, the mesh of the computational domain moves with the particle ¯uid velocity. In the Lagrangian description, the motion of the mesh may lead to an element entanglement; this description is preferred for problems with small motion. For problems involving free surface ¯ows, and problems with moving boundaries, it is necessary to have a formulation, which allows the mesh to follow the moving boundaries, and preserve the element shape. This has been achieved by the arbitrary Lagrangian±Eulerian (ALE) description. The ALE formulation has been introduced in the ®nite dierence literature by Noh et al. [1] and Hirt et al. [2]. Analogous ®nite element ALE formulation for incompressible viscous ¯ows has been introduced by Hughes et al. [3], Liu et al. [4] and Belytschko and Kennedy [5]. The freedom in moving the mesh oered by the ALE formulation is
*
Corresponding author. E-mail address: [email protected] (M. Souli).
0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 3 1 3 - 9
452
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
very attractive. The purpose of the ALE algorithm is to provide a capability for automatic reasoning that conserves the regularity of the computational mesh. Several methods have been used to solve free surface problems. One method, shape optimization, has been developed for ¯uid problems by Pironneau [6]. The method has been used for to instability problems in [7], and for water wave problems in [8]. Shape optimization method involves the dierentiation of a functional energy, with respect to the domain for continuous optimization, or with respect to the nodes of the domain for discrete optimization. For complex problems, ®nite dierence methods are used for the dierentiation, see [8]. This method is very successful for steady problems. For transient problems, the free surface function is governed by a parabolic type equation, see [9]. The solution of the equation gives the position of the free surface nodes. In order to move the internal nodes, a linear distribution of the nodes along the normal direction of the free surface is used. When using the ALE method, the solution of the free surface problem involves the solution of the ¯uid and mesh equations. For the ¯uid equations, several ®nite element formulations have been developed to solve the Navier±Stokes equations. For convection dominated ¯ows, the characteristic-Galerkin formulation was applied with success for incompressible ¯ows by Glowinski and Pironneau [10]. Another method, the streamline-upwind/Petrov±Galerkin (SUPG) method, introduced by Hughes and Brooks [11], has been applied ®rst to a linear scalar-diusion and to uncompressible Navier±Stokes equations. In this paper, we use the Galerkin/least-square (GLS) formulation to solve the ¯uid equations. The GLS method has been introduced by Hughes et al. [12], and has been successfully applied to the compressible Navier±Stokes equations by Shakib [13], and to incompressible ¯ows by Johnson et al. [14]. The outline of this paper is arranged as follows. In Section 2, a general theory of the ALE method and free surface description is introduced. We describe the partial dierential equations governing the mesh motions as well as the boundary conditions for the free surface, which evolves the normal velocity of the ¯uid. The ¯uid equations are coupled to the mesh equations through the free surface boundary conditions. A staggered method is used to solve the coupled problems. In Section 3, we describe the ®nite element formulation and time integration method for the ¯uid equations. Section 4 is devoted to the ®nite element formulation of the mesh equations and the free surface boundary conditions. This leads to a nonlinear system of dierential equations which is solved by a predictor±corrector algorithm. In order to conserve the mass when solving for the mesh motion, accurate approximations of the normal direction at the free surface are required. In Section 5, numerical results that indicate the eectiveness of the approach are presented. 2. The ALE and free surface description 2.1. The ALE description The ALE description for incompressible viscous ¯ows has been developed by Hughes et al. [3], to solve free surface ¯ows and ¯uid±structure interaction problems. A general kinematics theory was developed in [3], which serves as the basis of the Lagrangian±Eulerian description. For this purpose, the authors de®ne three domains in space, and mappings from one domain to the other. The ®rst one, called the spatial domain, is considered as the domain on which the ¯uid problem is posed. The spatial domain is generally in motion, because of moving boundaries. The second domain, called the material domain, is to be thought of as the domain occupied at time t 0 by the material particles which occupy the spatial domain at time t. The third domain, called the reference domain, is de®ned as a ®xed domain throughout. From these domain descriptions, we can see that the Eulerian description is obtained when the spatial domain coincides with the reference domain, whereas the Lagrangian reference is obtained when the material domain coincides with the reference domain. Both the material and spatial domains are generally in motion with respect to the reference domain; it is convenient to express the material time derivative of a physical property / in the reference con®guration /_ /;t c r/;
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
453
where /_ is the material time derivative, and /;t is the time derivative when freezing coordinates in the reference domain, c is the convective velocity cv
vmesh ;
where v is the ¯uid velocity, and vmesh is the mesh velocity. In the Eulerian description, the mesh velocity is zero, vmesh 0, whereas in the Lagrangian description vmesh v, and c 0. 2.2. The ALE and free surface description In the ALE formulation, the mesh nodes move with an arbitrary velocity. The choice of the mesh velocity constitutes one of the major problems with the ALE description. Dierent techniques have been developed for updating the mesh in a ¯uid motion, depending on the ¯uid domain. For problems de®ned in simple domains, the mesh velocity can be deduced through a uniform or nonuniform distribution of the nodes along straight lines ending at the moving boundaries. This technique has been used in [3,4,6], for water wave problems. For general computational domains, the mesh velocity is computed through partial dierential equations, with appropriate boundary conditions, see [15]. In [16], for a two-dimensional problem, Masud introduced the following equation for the mesh displacement ®eld d: div 1 srd 0; where s is a nondimensional function given on each element by s
Dmax Dmin Delem
Dmin represents the area of the current element. The parameter s is designed in order to prevent distortion of small elements. In this paper, a more general partial dierential equation than the Laplacian operator is used. A nonlinear elasticity model is introduced to simulate the mesh motion. The computational domain is considered as an elastic body and is deformed by external loads applied throughtout the moving boundary. The partial dierential equations governing the displacement ®eld of the mesh are given by div r 0
1
with appropriate boundary conditions, r is the Cauchy stress tensor. To encompass most of the ¯uid problems, it is necessary to consider two cases, the ®rst one concerns an incompressible domain, usually used for contained ¯uid problems. Since the mass of the ¯uid is conserved, the domain of the incompressible ¯uid should be conserved, when solving for the mesh motion. This situation arises in sloshing tank problems. The second case concerns problems in which the volume of the computational domain is not necessarily conserved even though the ¯uid equations satisfy mass conservation, and the total boundary ¯ux is zero. This situation arises for draining tank problem in which the mesh has to be compressed. For the ®rst case, in which the mesh is viewed as incompressible, we consider a linear elasticity model, where the constitutive relation is given by r C e:
2
The constitutive matrix C is element dependent. C is constructed such that smaller elements undergo lesser deformation. The constitutive matrix coef®cients are inversely proportional to the element size. The design of C helps to keep the mesh re®ned in the areas where small scale effects are of interest, as in boundary layer problems. The second case, the compressible formulation for the mesh motion, is given by a nonlinear neo-Hookean elasticity model. For this problem, the strain energy function is given by
454
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
1 w e k J 2
1 2 1 l I1 2
3
l log J ;
3
where k and l are the material parameters which are equivalent to the Lame coecients for small deformations, and I1 and J are the ®rst and third invariants of the deformation matrix. In order to solve the mesh Eq. (1), boundary conditions need to be applied: (a) Fixed boundary conditions. The mesh displacement is zero d mesh 0
on S0 :
4
(b) Moving boundary conditions. This is the case of a moving body in the ¯uid or moving walls. For this boundary, the mesh displacement is prescribed d mesh g
on Sg :
5
(c) Free surface boundary conditions. No ¯uid particle can cross the free surface. We express this condition by setting the normal mesh velocity at the free surface equal to the normal ¯uid velocity vmesh ~ n 0;
v
6
where v is the ¯uid velocity, and ~ n is the unit normal vector. The mesh equations and the appropriate boundary conditions are given by div r 0 d
mesh
0 on S0 fixed boundary;
d mesh g v
in X; on Sg moving boundary;
n0 vmesh ~
7
on S free surface:
Remark. In most free surface problems, the dynamic of the free surface g t; x1 ; x2 is governed by a hyperbolic type equation, see [9], described by og v rg 0; ot
8
where v is the ¯uid particle velocity. The numerical solution of problem (8) gives the location of the free surface nodes. 2.3. ALE description of the Navier±Stokes equations Let X 2 Rn , n 2 or 3, the domain occupied by the ¯uid particles, and let C denote its boundary. The equations of mass and momentum conservation for an incompressible Newtonian ¯uid in the ALE form are given by: div v 0; q
ov q v ot
9 vmesh rv div r f ;
10
where q is the density, and r is the total Cauchy stress given by r
T
p Id l rv rv ;
where p is the pressure and l is the dynamic viscosity. The boundary conditions are introduced as follows. The part of the boundary at which the velocity is assumed to be speci®ed is denoted by Cg . This is an in¯ow boundary condition vg
on Cg :
11
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
455
The traction boundary conditions associated with Eqs. (9) and (10) are the conditions on stress components. These conditions are assumed to be imposed on the remaining part of the boundary r ~ nh
on Ch :
12
The homogeneous counterpart of (12) corresponds to traction free conditions, and is often at the out¯ow boundaries. This condition is referred to as a natural boundary condition r ~ n0
on C0 :
13
As initial condition, a divergence free velocity ®eld u0 x is speci®ed over the domain at t 0 v x; 0 v0 x:
14
3. Finite element formulation for ¯uid equations 3.1. The space discretization The Navier±Stokes Eqs. (9) and (10) with appropriate boundary conditions (11)±(13), and initial condition (14) are discretized in space using a GLS formulation. In this formulation, a least-square form of the residual is added to the Galerkin equations. These terms enhance the stability of the Galerkin method without losing accuracy. In the GLS formulation we are using both velocity and pressure which are approximated with linear interpolation Q1 . This equal interpolation has been shown to be stable within the formulation, and circumvent the Babuska±Brezzi condition, see [17]. The ®nite element discretization of the domain X is achieved by dividing it into elements. Associated with the discretization, we de®ne the following ®nite element interpolation function spaces for the velocity and pressure: 3 V v 2 H 1 X ; v j Cg g; v j K 2 Q1 g; 3 W w 2 H 1 X ; w j Cg 0; v j K 2 Q1 g; 3 P p 2 H 1 X ; p j K 2 Q1 g: The ®nite element formulation is written as follows: Find u 2 V and p 2 P , such that 8w 2 W and 8q 2 P , we have Z I v rq dx v n q ds 0; Z X
X
C
q v;t v rvw dx
Z
I X
r rw dx
Z Ch
h w ds
X
f v dx
15 XZ K
K
I W s I V
f dv 0:
The last integral in (15) is the least-square term. I is de®ned as the incompressible Navier±Stokes operator I V q A~0;t V A~i V;i
K~ij;i V;i j :
This operator is applied to the vector V v; p=q and the weight function W w; q. The 4 4 matrices A~0 ; A~i ; K~i;j are given by 0 1 1 0 0 0 B0 1 0 0C 0 ei B C ~ ~ ~ A0 B C ; Ai v i A 0 @0 0 1 0A eti 0 0 0
0
1
16
456
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
and K~ij 0 for 0 l B0 K~ii B @0 0
i 6 j and 0 l 0 0
0 0 l 0
1 0 0C C: 0A l
Remark. 1. The least-square integral is only de®ned on the interior of the elements. 2. An integration by parts is used for the mass equation, in order to enforce mass conservation. This can be seen by considering w 0 and q 1 in Eq. (15), which gives a zero total mass ¯ux I u n ds 0: 17 C
3. The construction of the s matrix has been developed for compressible ¯ows by Shakib [13]. For the onedimensional advective diusive problem, an optimal construction of the parameter s is given in [13] by i1=2 h 2 2 18 s 2=Dt 2kuk2 =h 9 4m=h2 ; where m is the kinematic viscosity, Dt is the time step, and h is the mesh size. For incompressible ¯ows, the s matrix is a multidimensional generalization of (18). 3.2. Time integration The discretization of the ®nite element formulation (15) leads to the following nonlinear discrete equations: M
dU N U 0; dt
19
where M is the mass matrix, U the vector of nodal values, and N U is a nonlinear vector of U. To solve the semi-discrete system (19), we develop a predictor±corrector algorithm to reduce the nonlinear system to a sequence on linear systems. In order to introduce numerical dissipation without degrading the order of accuracy, Hilbert, Hughes and Taylor developed the a-method [18], referred to as HHT method. In the HHT method, the ®nite dierence formulas of the Newmark method are retained, whereas the semi-discrete equation is modi®ed as follows: M U_ n1 N aUn1 1
aUn Fna ;
20
where Fna aFn1 1
aFn :
The following equation completes the method: Un1 Un Dt 1
cU_ cU_ n1 ;
21
where c 2 0; 1 is a relaxation parameter. The parameters a and c are selected to obtain stability and accuracy [19]. A nonlinear algebraic problem may be obtained for U_ n1 by eliminating Un1 from Eqs. (20) and (21), The nonlinear algebraic problem is solved using a predictor±corrector algorithm. There are several possible implementations of the algorithm. For ¯uid problems, we employ the so-called velocity form. The ®nal solution of the velocity at the previous time step is used as predictor value for the current time step. The predictors for U and U_ are given by
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466 0 Un1 Un ; 0 _ Un1 1 1=cU_ n :
457
22
The system is solved using Newton's method starting with the predictors. The dierentiation of Eqs. (20) and (21) with respect to U_ n1 leads at each Newton's iteration to the following linear system: M acDtN~ Dain1
i M U_ n1
i N Una Fna ;
23
where i i Una aUn1 1
aUni
i and N~ oN =oU Una is the Jacobian matrix. Once the linear system (23) is solved for the incremental i acceleration Dan1 , the corrector phase is given by i1 i1 U_ n1 Un1 Dain1 ; i1 i Un1 Un1 cDtain1 :
24
Remark. 1. If the parameters a and c are selected such that a 2 2=3; 1 and c 1=2 3 2a, an unconditionally stable second-order accurate results, for more details, see [19]. 2. If a 1, the a-method is reduced to a Newmark method. In the HHT method decreasing a increases the amount of numerical dissipation. This can be exploited in quickly bringing about a transient solution to equilibrium. For steady solutions, one uses the parameters a 2=3 and c 5=6, which give an unconditionally stable and second-order algorithm. 3. The time increment Dt is computed based on the Courant±Friedricks±Levy (CFL) condition. 4. The convergence of the algorithms (20) and (21) is very sensitive to the choice of the Jacobian matrix N~ . It has been shown in [20], through numerical tests, that a consistent Jacobian of the left-hand side of (23) yields in general the fastest convergence to the steady-state solution, provided the initial guess is within the radius of convergence of Newton's algorithm. For most ¯uid problems, the use of the consistent Jacobian leads to poor stability in the initial phase of the time marching convergence. In order to improve stability, the advection and the least-square terms are assumed as constant in the process of dierentiating Eq. (20). 4. Finite element formulation for ALE equations 4.1. Space discretization In the absence of the free surface boundary condition, the boundary value problem can be studied as a minimization of the potential p. Let us de®ne the following spaces: o n 3 V d 2 H 1 X ; d j Sg d; d j S0 0 ; o n 3 W w 2 H 1 X ; w j Sg 0; w j S0 0 : Problem (7) is equivalent to the following problem. Find d 2 V , which minimizes the following strain energy: Z p d w g dx; X
where g d is the strain tensor for the displacement d, and w g the strain energy function.
25
458
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
The Dirichlet boundary conditions are embedded in the admissible space functions. From a numerical point of view, these conditions are usually solved by eliminating the equations related to these nodes. The free surface condition is a nonlinear condition, and can be satis®ed through penalization. To satisfy, we use an augmented Lagrangian formulation. In the augmented Lagrangian formulation, the Lagrange multiplier k is an additional unknown that we need to determine in addition to the ®eld variables. Let us introduce the new energy p1 I e p1 d p d kF vmesh F vmesh 2 dx: 26 2 S Throughout this paper, we will denote by d, u and v the mesh displacement, mesh velocity and ¯uid particle velocity, respectively, and by S the free surface. The function F u u v ~ n de®nes the constraint on the free surface S, e a penalty parameter and p is the potential corresponding to the elasticity problem (25). The variational equation emanating from Eq. (26) is given by Z I ow g d g dx kF;u e F F;u du dx 0; 27 og X S where dg
1 T rdw rdw 2
and
F;u
oF u : ou
We discretize the variational problem (27) using linear interpolation for the displacement d and nodal Dirac function for du. The ®nite element discretization of (27) leads to the following nonlinear system at each time step: K d
Ns X
kj eF uj Nj 0;
28
j1
Ns is the number of nodes on the free surface S, kj and uj are the nodal values of the Lagrangian parameters and the mesh velocity at the node j. The vectors Nj 2 RNdof , where Ndof is the total number of degrees of freedom of the system. These vectors are constructed in order to add the penalization terms to the appropriate equations in the linear system. For each node j, these vectors are de®ned by Nj 0; . . . ; 0; n1j ; n2j ; n3j ; 0; . . . ; 0 where nkj are the components of the unit normal vector nj at the node j. In taking the derivation of the constraint F u with respect to the mesh velocity u, we assume that the normal vector n is constant over the time step. 4.2. Free surface discretization The nonlinear system necessitates the values of the unit normal vector at the nodes on the free surface. The nodal value of the unit normal is not well de®ned when using linear interpolation, dierent approaches for de®ning the normal vector can be used. Mass conservation in the Navier±Stokes problem is very sensitive to the particular approach used for the normal. For this reason, it constitutes one of the major problems in the ALE description. In order to determine the best approach, we consider a situation, which arises in sloshing tank problems. In these problems, ¯uid mass can only pass through the free surface boundary, and thus the mass conservation is expressed by I v n ds 0: 29 S
For an incompressible ¯uid, the conservation of mass is equivalent to the conservation of the ¯uid volume. To conserve the volume of the ¯uid domain, the mesh velocity u should satisfy
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
459
Z S
u n ds 0:
30
The free surface condition is satis®ed on each node on the free surface through the augmented Lagrangian formulation uj
vj nj 0;
j 1; . . . ; Ns :
The normal vector nj at the node j is de®ned by I nj wj n ds;
31
32
wj is the shape function at the node j. From Eq. (32), the normal vector is well de®ned. Substituting (32) into (31) and summing over the free surface nodes, we get I v u n ds 0; 33 S
PNs PNs where u j1 uj wj and v j1 vj wj . Volume conservation is deduced from equation of mass and Eq. (33). 4.3. Time integration In order to obtain a semi-discrete system governing the mesh motion, with stability properties, let us introduce a mass matrix M and a damping matrix C de®ned by ~ M mK;
~ C cK;
34
where the parameters m and c are positive and will be determined later. K~ oK d=od is the tangent stiness. Adding the mass and the damping terms to the nonlinear system (28), we get at each time step the following system: Man1 Cun1 K dn1
Ns X kn1 eF un1 Nj 0;
35
j1
where an1 is the mesh acceleration and un1 is the mesh velocity. To keep down the proliferation of notations, we dropped the index j in the penalization terms kn1 and un1 . The nonlinear system (35) is completed by the following ®nite dierence formulae: Dt2 1 2ban 2b an1 ; 2 can c an1 :
dn1 dn Dt un un1 un Dt 1
36
The parameters a and b are selected for numerical stability. To solve the systems (35) and (36), we use the HHT method for numerical dissipation Man1 Cuna K dna
Ns X kna eF una Nj 0; j1
where
37
460
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
una a un1 1
aun ;
dna a dn1 1
adn ;
kna a kn1 1
akn :
In order to solve the nonlinear system (37), a predictor±corrector algorithm is used. For the mesh motion, the so-called displacement formulation is implemented. The ®nal displacement of the previous time step is used as a predictor value for the current time step. The predictors for the displacement, velocity and acceleration are given by 0 dn ; dn1
u0n1 1 a0n1
c=bun Dt 1 c=2ban ; 1 1 un 1 an : b Dt 2b
38
A nonlinear algebraic problem is obtained for the acceleration an1 by eliminating dn1 and un1 from Eqs. (37) and (38). At each time step, the nonlinear algebraic problem for an1 is solved using Newton's method starting from the predictor solution. In order to decrease the nonlinear eects of the system, we assume the normal vector constant over the time step. The normal vector is not updated for the nonlinear Newton's method. Using the expressions of the mass matrix M and the damping matrix C, we get the following linear system to solve: i Given ain1 ; uin1 ; dn1 and kn . First, we update the Lagrange multiplier, using Uzawa's method kin1 kn
eF uin1 :
Second, we solve the following linear system for the corrector phase: ( ) Ns X m cabDt acDt2 K~ acDte Nj Nj Dai1 Mai Ri j1
n1
n1
na ;
39
where i Rina Cuina K dna
Ns X j1
kina F uina Nj :
40
Once the linear system (39) is solved, the corrector phase is given by i i1 ai1 n1 an1 Dan1 ;
i1 un1 uin1 cDt ai1 n1 ;
i1 i dn1 dn1 bDt ai1 n1 :
The parameters m and c in (35) are selected such that the linear homogeneous part of the system (36) governing the mesh motion, given by ~ K ma acu ad 0
41
is critically damped. The linear system (41) can be considered as single degree of freedom equation m
d2 x dx ca ax 0: 2 dt dt
42
The solution of the single Eq. (42) is critically damped if ma2 c2
4ma 0
43
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
and is given by x t x 0 exp
ac at : 2m
461
44
The second relationship between the parameters m and c is deduced by assuming that the solution (44) decreases by 90% at the ®rst step, which gives acDt exp 0:1: 45 2m From relations (44) and (45), we get the values of the parameters m and c m
aDt2 log 0:1
2
;
c 2 m=a
1=2
:
Remark. If the parameters a; b and c are chosen such that a 2 2=3; 1, b 1=4 2 an unconditionally stable second-order accurate scheme results, see [19].
a2 and c 1=2 3
a,
5. Numerical examples To determine the eectiveness of the ALE method, we apply the formulation to three dierent engineering problems. The ®rst is a ®nite element problem of wave propagation. The numerical simulation of this problem has been carried out, ®rst by Hughes et al. [3], where the results are compared to the experimental ones, see [12]. The second example involves a large amplitude sloshing tank. Numerical results are compared to those given by Liu et al. [14]. In the last example, a rotation of ¯uid in a cylindrical container about its central axis is considered. Analytical solution exists for the steady problem, which consists of a rigid body rotation. 5.1. The solitary wave problem The solitary wave traveling in a rectangular channel is de®ned as a single elevation above a surrounding undisturbed water level. The wave is neither followed nor preceded by another elevation of the surface. The wave travels without change of shape and essentially constant velocity. Shallow water theory is used for the analytical approach. This problem has been studied by Hughes et al. [3]. The wave is generated by prescribing the velocity time history of the ¯uid along the left-hand boundary of the domain. The physical wave parameters are selected in keeping with the numerical simulation in [3], and the experimental wave generation facility in Keck Laboratory at the California Institute of Technology. The analytical time history of displacement of the wave is given by cjt i gh d 1 tanh 4 ; j d c g D p1=2 ; 1=2 3g ; j D where g 1, p 0:86, L 9949:094 the length of the channel, and D 10 the depth of the channel. The time increment Dt h=c is selected such that cfl 1, h is the element size. The wave-generation function is designed to produce a solitary wave, one in which nonlinear dispersive effects are balanced so that the wave propagated without distortion. The computational domain consists of 160 elements mesh in the horizontal direction, and two elements through the depth. As indicated in [3], by using two elements through the depth, most of the dispersion emanating from the tail of the wave can be avoided. The numerical results are in good
462
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
agreement with those given in [3]. The variation of the wave volume during the 150 time steps is less than 2%, and the amplitude decay is less than 3%. This is due to the fact that the GLS formulation models accurately the ¯uid incompressibility. The fact is that we use several passes for both the ¯uid and the mesh stagger, instead of one, with the mass and volume conservation being accurately satis®ed. This clearly shows that the numerical method is not aected by any arti®cial damping eect. In order to get less amplitude decay and conservation of the wave volume, Newmark parameters for the time integration are considered, thus the parameters a 1, b 1=2, c 1=2 are used for both ¯uid and ALE equations. Numerical experience shows that the time derivative in the least-square terms in Eq. (16) should be considered for transient problems; when this term is not included, a greater decay of wave amplitude is obtained. For steady-state problems, this term does not need to be included in the ®nite element formulation. The successful application of the present approach to this problem indicates the possibility of employing this technique to a wide variety of water wave problems where analytical methods are dicult, and sometimes impossible to derive. Fig. 1 shows the wave shape after 50, 100, 150 time steps. 5.2. The sloshing tank To show the adaptivity of the ALE formulation for large mesh distortion, the analysis of large amplitude sloshing in a rectangular tank has been carried out. The problem has been studied by Huerta and Liu [4]. A schematic view of the tank and the boundary conditions are shown in Fig. 2. The width of the tank is L 0:8 m and the depth D 0:3 m. The tank is excited by a sinusoidal acceleration a A g sin xt;
46
where g; x and t are the gravity, the circular frequency and time, respectively. A is an arbitrary constant 1=2 governing the amplitude of the acceleration. The Reynolds number based on the reference velocity g D 5 is 5 10 . Gravity acts downwards and the amplitude of the motion is A 0:001. Using these parameters, experimental results show that the resonant frequency of the tank is x 0:88 Hz. The time increment Dt p=30x, and the sloshing period is T 2p=x. The ®nite element mesh consists of 441 21 21 elements. In order to conserve mass during the mesh motion, one uses the following strategy at each time step:
Fig. 1.
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
463
Fig. 2.
First, we run the ¯uid problem to convergence, which insures mass conservation I v n ds 0; S
47
Second, we run the mesh problem to convergence, and from free surface boundary conditions and mass conservation, we have I u n ds 0: 48 S
Eq. (48) indicates that at each time step, the mesh motion conserves the computational volume and thus mass conservation is satis®ed. Fig. 3 shows the vertical displacement of the nodes at the wall tank, as a function of time. The results are in good agreement with the results given in [4,21].
Fig. 3.
464
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
Fig. 4.
Fig. 5.
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
465
5.3. The rotating bucket Let us consider a ¯uid in a cylindrical container. The cylinder is rotating about its central axis, a schematic view is shown in Fig. 4. We assume that the container has been rotating at constant velocity long enough for the ¯uid to have attained rigid rotation. In the rotating frame, the ¯uid has no-slip boundary conditions at the wall and at the bottom of the cylinder. The ¯uid is subjected to gravity body force and the external force due to the centripetal acceleration f q x x xxr, where r is the distance to the central axis, and x is the rotational velocity. At steady-state, the free surface de®ned by the free surface condition (6) is a parabolic. Fig. 5 shows the graph of exact and numerical free surfaces. This example illustrates that the numerical integration parameters a 2=3, b 4=9 and c 9=4 lead to an accurate solution for the steady-state problem.
6. Conclusion In this paper, an ALE free surface formulation for viscous incompressible ¯ow has been presented. The GLS method is used for the ®nite element spatial discretization of the Navier±Stokes equations. A predictor±corrector scheme with numerical damping is adopted for time integration. For the mesh motion, elasticity type equations, with appropriate boundary conditions are used. Two types of constitutive models are used for the mesh motion. The ®rst one is a neo-Hookean model for the incompressible ¯ows, and draining tank problem. For the mesh motion, an appropriate discretization of the normal vector is used, in order to satisfy mass conservation in the computational domain. Numerical examples show that the ALE method is appropriate for steady and transient free surface problems.
References [1] W.F. Noh, A time dependent two space dimensional coupled Eulerian Lagrangian code, Meth. Comput. Phys. 3 (1964). [2] C.W. Hirt, A.A. Amsden, H.K. Cook, An arbitrary Lagrangian±Eulerian method for all ¯ow speeds, J. Comput. Phys. 14 (1974) 227±253. [3] T.J.R. Hughes, W.K. Liu, T.K. Zimmerman, Lagrangian±Eulerian ®nite element formulation for viscous ¯ows, Comput. Methods Appl. Mech. Engrg. 29 (1981) 329±349. [4] A. Huerta, W.K. Liu, Viscous ¯ows for large free surface motion, Comput. Methods Appl. Mech. Engrg. 69 (1988) 277±324. [5] T. Belytschko, D.F. Flanagan, J.M. Kennedy, Finite element method with user-controlled meshes for ¯uid±structure interactions, Comput. Methods Appl. Mech. Engrg. 33 (1982) 689±723. [6] O. Pironneau, Shape Optimization for Elliptic Problems, Springer, Berlin, 1983. [7] M. Souli, Y. Demay, A. Habbal, Finite element study for the draw resonance instability, Eur. J. Mech./Fluids 11 (6) (1992). [8] M. Souli, J.P. Zolesio, Shape derivative of discretized problems, Comput. Methods Appl. Mech. Engrg. 108 (1993). [9] B. Ramaswany, M. Kawahara, Arbitrary Lagrangian±Eulerian ®nite element method for unsteady convective incompressible viscous free surface ¯uid ¯ow, Int. J. Numer. Methods Fluids 7 (1987) 1053±1075. [10] R. Glowinski, O. Pironneau, Finite element methods for Navier±Stokes equations, Ann. Rev. Fluid Mech. 24 (1992) 167±204. [11] A.N. Brooks, T.J.R. Hughes, Streamline upwinded/Petrov±Galerkin methods for advection dominated ¯ows with particular emphasis on the incompressible Navier±Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199±259. [12] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new ®nite element formulation for computational ¯uid dynamics, Comput. Methods Appl. Mech. Engrg. 54 (1986) 223±234. [13] F. Shakib, Finite element analysis of the compressible Euler and Navier±Stokes equations, Ph.D. Thesis, Stanford University, 1988. [14] C. Johnson, J. Saranem, Diusion methods for the incompressible Euler and Navier±Stokes equations, Math. Comput. 47 (1986) 1±18. [15] T. Tezduyar, J. Liou, M. Behr, New strategy for FEM computations involving moving boundaries and interfaces, I, USMI Report 90/69. [16] A. Masud, A space-time ®nite element method for ¯uid±structure interactions, Ph.D. Thesis, Stanford University, 1993. [17] T.J.R. Hughes, L.P. Franca, M. Balestra, A new ®nite element formulation for computational ¯uid dynamics. V. Circumventing the Babuska±Brezzi condition, a stable Petrov±Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85±99.
466
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
[18] H.M. Hilbert, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms for structural dynamics, Engrg. Struct. Dyn. 5 (1977) 283±292. [19] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewoods Clis, NJ, 1987. [20] Z. Johan, T.J.R. Hughes, F. Shakib, A matrix free implicit iterative solver for compressible ¯ow problems. [21] D. Gorin, Ph.D., California Institute of Technology, 1979.
www.elsevier.com/locate/cma
Arbitrary Lagrangian±Eulerian and free surface methods in ¯uid mechanics M. Souli a
a,*
, J.P. Zolesio
b
Lab. de Mechanique, Universit e des Sciences et Technologie de Lille 1, Bd Paul Langevin, Cite Scienti®que, 59655 Lille, France b Ecole des Mines de Paris, Sophia Antipolis, Route des Lucioles, Valbonne, France Received 25 March 2000; received in revised form 25 August 2000; accepted 30 May 2001
Abstract A computational procedure is developed to solve problems of viscous incompressible ¯ows under large free surface motions. The arbitrary Lagrangian±Eulerian (ALE) method is used to move the free surface nodes as well as the internal nodes. The coupling of the mesh motion equations and the ¯uid equations is essentially done through the free surface boundary conditions. An elastic model with free surface boundary conditions is used to model the mesh motion. The boundary condition at the free surface is deduced from the fact that no ¯uid particle can cross the free surface. The ¯uid problem is solved using a Galerkin/least-square (GLS) ®nite element formulation. Ó 2001 Elsevier Science B.V. All rights reserved.
1. Introduction In continuum mechanics, two descriptions are considered for the motion in a continuum media: The ®rst is the Eulerian description, where we focus attention on a particular volume in space. The volume is ®xed with respect to a laboratory frame, and we study the ¯uid as it passes through the ®xed volume. The description is one in which the ¯uid is continuously renewed inside the volume, the Eulerian description is not the simplest in which the basic equations of ¯uid motion can be formulated. A convective term is introduced to express the material time derivative in the reference con®guration. The convective term gives a nonsymmetrical form of the Galerkin formulation. Since the computational domain is ®xed, the Eulerian description has the advantage of preserving the mesh regularity. The second is the Lagrangian description, in which we identify and follow a particular region of ¯uid. The volume of ¯uid changes in shape, while the total mass remains constant. In the Lagrangian description, the mesh of the computational domain moves with the particle ¯uid velocity. In the Lagrangian description, the motion of the mesh may lead to an element entanglement; this description is preferred for problems with small motion. For problems involving free surface ¯ows, and problems with moving boundaries, it is necessary to have a formulation, which allows the mesh to follow the moving boundaries, and preserve the element shape. This has been achieved by the arbitrary Lagrangian±Eulerian (ALE) description. The ALE formulation has been introduced in the ®nite dierence literature by Noh et al. [1] and Hirt et al. [2]. Analogous ®nite element ALE formulation for incompressible viscous ¯ows has been introduced by Hughes et al. [3], Liu et al. [4] and Belytschko and Kennedy [5]. The freedom in moving the mesh oered by the ALE formulation is
*
Corresponding author. E-mail address: [email protected] (M. Souli).
0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 3 1 3 - 9
452
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
very attractive. The purpose of the ALE algorithm is to provide a capability for automatic reasoning that conserves the regularity of the computational mesh. Several methods have been used to solve free surface problems. One method, shape optimization, has been developed for ¯uid problems by Pironneau [6]. The method has been used for to instability problems in [7], and for water wave problems in [8]. Shape optimization method involves the dierentiation of a functional energy, with respect to the domain for continuous optimization, or with respect to the nodes of the domain for discrete optimization. For complex problems, ®nite dierence methods are used for the dierentiation, see [8]. This method is very successful for steady problems. For transient problems, the free surface function is governed by a parabolic type equation, see [9]. The solution of the equation gives the position of the free surface nodes. In order to move the internal nodes, a linear distribution of the nodes along the normal direction of the free surface is used. When using the ALE method, the solution of the free surface problem involves the solution of the ¯uid and mesh equations. For the ¯uid equations, several ®nite element formulations have been developed to solve the Navier±Stokes equations. For convection dominated ¯ows, the characteristic-Galerkin formulation was applied with success for incompressible ¯ows by Glowinski and Pironneau [10]. Another method, the streamline-upwind/Petrov±Galerkin (SUPG) method, introduced by Hughes and Brooks [11], has been applied ®rst to a linear scalar-diusion and to uncompressible Navier±Stokes equations. In this paper, we use the Galerkin/least-square (GLS) formulation to solve the ¯uid equations. The GLS method has been introduced by Hughes et al. [12], and has been successfully applied to the compressible Navier±Stokes equations by Shakib [13], and to incompressible ¯ows by Johnson et al. [14]. The outline of this paper is arranged as follows. In Section 2, a general theory of the ALE method and free surface description is introduced. We describe the partial dierential equations governing the mesh motions as well as the boundary conditions for the free surface, which evolves the normal velocity of the ¯uid. The ¯uid equations are coupled to the mesh equations through the free surface boundary conditions. A staggered method is used to solve the coupled problems. In Section 3, we describe the ®nite element formulation and time integration method for the ¯uid equations. Section 4 is devoted to the ®nite element formulation of the mesh equations and the free surface boundary conditions. This leads to a nonlinear system of dierential equations which is solved by a predictor±corrector algorithm. In order to conserve the mass when solving for the mesh motion, accurate approximations of the normal direction at the free surface are required. In Section 5, numerical results that indicate the eectiveness of the approach are presented. 2. The ALE and free surface description 2.1. The ALE description The ALE description for incompressible viscous ¯ows has been developed by Hughes et al. [3], to solve free surface ¯ows and ¯uid±structure interaction problems. A general kinematics theory was developed in [3], which serves as the basis of the Lagrangian±Eulerian description. For this purpose, the authors de®ne three domains in space, and mappings from one domain to the other. The ®rst one, called the spatial domain, is considered as the domain on which the ¯uid problem is posed. The spatial domain is generally in motion, because of moving boundaries. The second domain, called the material domain, is to be thought of as the domain occupied at time t 0 by the material particles which occupy the spatial domain at time t. The third domain, called the reference domain, is de®ned as a ®xed domain throughout. From these domain descriptions, we can see that the Eulerian description is obtained when the spatial domain coincides with the reference domain, whereas the Lagrangian reference is obtained when the material domain coincides with the reference domain. Both the material and spatial domains are generally in motion with respect to the reference domain; it is convenient to express the material time derivative of a physical property / in the reference con®guration /_ /;t c r/;
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
453
where /_ is the material time derivative, and /;t is the time derivative when freezing coordinates in the reference domain, c is the convective velocity cv
vmesh ;
where v is the ¯uid velocity, and vmesh is the mesh velocity. In the Eulerian description, the mesh velocity is zero, vmesh 0, whereas in the Lagrangian description vmesh v, and c 0. 2.2. The ALE and free surface description In the ALE formulation, the mesh nodes move with an arbitrary velocity. The choice of the mesh velocity constitutes one of the major problems with the ALE description. Dierent techniques have been developed for updating the mesh in a ¯uid motion, depending on the ¯uid domain. For problems de®ned in simple domains, the mesh velocity can be deduced through a uniform or nonuniform distribution of the nodes along straight lines ending at the moving boundaries. This technique has been used in [3,4,6], for water wave problems. For general computational domains, the mesh velocity is computed through partial dierential equations, with appropriate boundary conditions, see [15]. In [16], for a two-dimensional problem, Masud introduced the following equation for the mesh displacement ®eld d: div 1 srd 0; where s is a nondimensional function given on each element by s
Dmax Dmin Delem
Dmin represents the area of the current element. The parameter s is designed in order to prevent distortion of small elements. In this paper, a more general partial dierential equation than the Laplacian operator is used. A nonlinear elasticity model is introduced to simulate the mesh motion. The computational domain is considered as an elastic body and is deformed by external loads applied throughtout the moving boundary. The partial dierential equations governing the displacement ®eld of the mesh are given by div r 0
1
with appropriate boundary conditions, r is the Cauchy stress tensor. To encompass most of the ¯uid problems, it is necessary to consider two cases, the ®rst one concerns an incompressible domain, usually used for contained ¯uid problems. Since the mass of the ¯uid is conserved, the domain of the incompressible ¯uid should be conserved, when solving for the mesh motion. This situation arises in sloshing tank problems. The second case concerns problems in which the volume of the computational domain is not necessarily conserved even though the ¯uid equations satisfy mass conservation, and the total boundary ¯ux is zero. This situation arises for draining tank problem in which the mesh has to be compressed. For the ®rst case, in which the mesh is viewed as incompressible, we consider a linear elasticity model, where the constitutive relation is given by r C e:
2
The constitutive matrix C is element dependent. C is constructed such that smaller elements undergo lesser deformation. The constitutive matrix coef®cients are inversely proportional to the element size. The design of C helps to keep the mesh re®ned in the areas where small scale effects are of interest, as in boundary layer problems. The second case, the compressible formulation for the mesh motion, is given by a nonlinear neo-Hookean elasticity model. For this problem, the strain energy function is given by
454
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
1 w e k J 2
1 2 1 l I1 2
3
l log J ;
3
where k and l are the material parameters which are equivalent to the Lame coecients for small deformations, and I1 and J are the ®rst and third invariants of the deformation matrix. In order to solve the mesh Eq. (1), boundary conditions need to be applied: (a) Fixed boundary conditions. The mesh displacement is zero d mesh 0
on S0 :
4
(b) Moving boundary conditions. This is the case of a moving body in the ¯uid or moving walls. For this boundary, the mesh displacement is prescribed d mesh g
on Sg :
5
(c) Free surface boundary conditions. No ¯uid particle can cross the free surface. We express this condition by setting the normal mesh velocity at the free surface equal to the normal ¯uid velocity vmesh ~ n 0;
v
6
where v is the ¯uid velocity, and ~ n is the unit normal vector. The mesh equations and the appropriate boundary conditions are given by div r 0 d
mesh
0 on S0 fixed boundary;
d mesh g v
in X; on Sg moving boundary;
n0 vmesh ~
7
on S free surface:
Remark. In most free surface problems, the dynamic of the free surface g t; x1 ; x2 is governed by a hyperbolic type equation, see [9], described by og v rg 0; ot
8
where v is the ¯uid particle velocity. The numerical solution of problem (8) gives the location of the free surface nodes. 2.3. ALE description of the Navier±Stokes equations Let X 2 Rn , n 2 or 3, the domain occupied by the ¯uid particles, and let C denote its boundary. The equations of mass and momentum conservation for an incompressible Newtonian ¯uid in the ALE form are given by: div v 0; q
ov q v ot
9 vmesh rv div r f ;
10
where q is the density, and r is the total Cauchy stress given by r
T
p Id l rv rv ;
where p is the pressure and l is the dynamic viscosity. The boundary conditions are introduced as follows. The part of the boundary at which the velocity is assumed to be speci®ed is denoted by Cg . This is an in¯ow boundary condition vg
on Cg :
11
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
455
The traction boundary conditions associated with Eqs. (9) and (10) are the conditions on stress components. These conditions are assumed to be imposed on the remaining part of the boundary r ~ nh
on Ch :
12
The homogeneous counterpart of (12) corresponds to traction free conditions, and is often at the out¯ow boundaries. This condition is referred to as a natural boundary condition r ~ n0
on C0 :
13
As initial condition, a divergence free velocity ®eld u0 x is speci®ed over the domain at t 0 v x; 0 v0 x:
14
3. Finite element formulation for ¯uid equations 3.1. The space discretization The Navier±Stokes Eqs. (9) and (10) with appropriate boundary conditions (11)±(13), and initial condition (14) are discretized in space using a GLS formulation. In this formulation, a least-square form of the residual is added to the Galerkin equations. These terms enhance the stability of the Galerkin method without losing accuracy. In the GLS formulation we are using both velocity and pressure which are approximated with linear interpolation Q1 . This equal interpolation has been shown to be stable within the formulation, and circumvent the Babuska±Brezzi condition, see [17]. The ®nite element discretization of the domain X is achieved by dividing it into elements. Associated with the discretization, we de®ne the following ®nite element interpolation function spaces for the velocity and pressure: 3 V v 2 H 1 X ; v j Cg g; v j K 2 Q1 g; 3 W w 2 H 1 X ; w j Cg 0; v j K 2 Q1 g; 3 P p 2 H 1 X ; p j K 2 Q1 g: The ®nite element formulation is written as follows: Find u 2 V and p 2 P , such that 8w 2 W and 8q 2 P , we have Z I v rq dx v n q ds 0; Z X
X
C
q v;t v rvw dx
Z
I X
r rw dx
Z Ch
h w ds
X
f v dx
15 XZ K
K
I W s I V
f dv 0:
The last integral in (15) is the least-square term. I is de®ned as the incompressible Navier±Stokes operator I V q A~0;t V A~i V;i
K~ij;i V;i j :
This operator is applied to the vector V v; p=q and the weight function W w; q. The 4 4 matrices A~0 ; A~i ; K~i;j are given by 0 1 1 0 0 0 B0 1 0 0C 0 ei B C ~ ~ ~ A0 B C ; Ai v i A 0 @0 0 1 0A eti 0 0 0
0
1
16
456
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
and K~ij 0 for 0 l B0 K~ii B @0 0
i 6 j and 0 l 0 0
0 0 l 0
1 0 0C C: 0A l
Remark. 1. The least-square integral is only de®ned on the interior of the elements. 2. An integration by parts is used for the mass equation, in order to enforce mass conservation. This can be seen by considering w 0 and q 1 in Eq. (15), which gives a zero total mass ¯ux I u n ds 0: 17 C
3. The construction of the s matrix has been developed for compressible ¯ows by Shakib [13]. For the onedimensional advective diusive problem, an optimal construction of the parameter s is given in [13] by i1=2 h 2 2 18 s 2=Dt 2kuk2 =h 9 4m=h2 ; where m is the kinematic viscosity, Dt is the time step, and h is the mesh size. For incompressible ¯ows, the s matrix is a multidimensional generalization of (18). 3.2. Time integration The discretization of the ®nite element formulation (15) leads to the following nonlinear discrete equations: M
dU N U 0; dt
19
where M is the mass matrix, U the vector of nodal values, and N U is a nonlinear vector of U. To solve the semi-discrete system (19), we develop a predictor±corrector algorithm to reduce the nonlinear system to a sequence on linear systems. In order to introduce numerical dissipation without degrading the order of accuracy, Hilbert, Hughes and Taylor developed the a-method [18], referred to as HHT method. In the HHT method, the ®nite dierence formulas of the Newmark method are retained, whereas the semi-discrete equation is modi®ed as follows: M U_ n1 N aUn1 1
aUn Fna ;
20
where Fna aFn1 1
aFn :
The following equation completes the method: Un1 Un Dt 1
cU_ cU_ n1 ;
21
where c 2 0; 1 is a relaxation parameter. The parameters a and c are selected to obtain stability and accuracy [19]. A nonlinear algebraic problem may be obtained for U_ n1 by eliminating Un1 from Eqs. (20) and (21), The nonlinear algebraic problem is solved using a predictor±corrector algorithm. There are several possible implementations of the algorithm. For ¯uid problems, we employ the so-called velocity form. The ®nal solution of the velocity at the previous time step is used as predictor value for the current time step. The predictors for U and U_ are given by
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466 0 Un1 Un ; 0 _ Un1 1 1=cU_ n :
457
22
The system is solved using Newton's method starting with the predictors. The dierentiation of Eqs. (20) and (21) with respect to U_ n1 leads at each Newton's iteration to the following linear system: M acDtN~ Dain1
i M U_ n1
i N Una Fna ;
23
where i i Una aUn1 1
aUni
i and N~ oN =oU Una is the Jacobian matrix. Once the linear system (23) is solved for the incremental i acceleration Dan1 , the corrector phase is given by i1 i1 U_ n1 Un1 Dain1 ; i1 i Un1 Un1 cDtain1 :
24
Remark. 1. If the parameters a and c are selected such that a 2 2=3; 1 and c 1=2 3 2a, an unconditionally stable second-order accurate results, for more details, see [19]. 2. If a 1, the a-method is reduced to a Newmark method. In the HHT method decreasing a increases the amount of numerical dissipation. This can be exploited in quickly bringing about a transient solution to equilibrium. For steady solutions, one uses the parameters a 2=3 and c 5=6, which give an unconditionally stable and second-order algorithm. 3. The time increment Dt is computed based on the Courant±Friedricks±Levy (CFL) condition. 4. The convergence of the algorithms (20) and (21) is very sensitive to the choice of the Jacobian matrix N~ . It has been shown in [20], through numerical tests, that a consistent Jacobian of the left-hand side of (23) yields in general the fastest convergence to the steady-state solution, provided the initial guess is within the radius of convergence of Newton's algorithm. For most ¯uid problems, the use of the consistent Jacobian leads to poor stability in the initial phase of the time marching convergence. In order to improve stability, the advection and the least-square terms are assumed as constant in the process of dierentiating Eq. (20). 4. Finite element formulation for ALE equations 4.1. Space discretization In the absence of the free surface boundary condition, the boundary value problem can be studied as a minimization of the potential p. Let us de®ne the following spaces: o n 3 V d 2 H 1 X ; d j Sg d; d j S0 0 ; o n 3 W w 2 H 1 X ; w j Sg 0; w j S0 0 : Problem (7) is equivalent to the following problem. Find d 2 V , which minimizes the following strain energy: Z p d w g dx; X
where g d is the strain tensor for the displacement d, and w g the strain energy function.
25
458
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
The Dirichlet boundary conditions are embedded in the admissible space functions. From a numerical point of view, these conditions are usually solved by eliminating the equations related to these nodes. The free surface condition is a nonlinear condition, and can be satis®ed through penalization. To satisfy, we use an augmented Lagrangian formulation. In the augmented Lagrangian formulation, the Lagrange multiplier k is an additional unknown that we need to determine in addition to the ®eld variables. Let us introduce the new energy p1 I e p1 d p d kF vmesh F vmesh 2 dx: 26 2 S Throughout this paper, we will denote by d, u and v the mesh displacement, mesh velocity and ¯uid particle velocity, respectively, and by S the free surface. The function F u u v ~ n de®nes the constraint on the free surface S, e a penalty parameter and p is the potential corresponding to the elasticity problem (25). The variational equation emanating from Eq. (26) is given by Z I ow g d g dx kF;u e F F;u du dx 0; 27 og X S where dg
1 T rdw rdw 2
and
F;u
oF u : ou
We discretize the variational problem (27) using linear interpolation for the displacement d and nodal Dirac function for du. The ®nite element discretization of (27) leads to the following nonlinear system at each time step: K d
Ns X
kj eF uj Nj 0;
28
j1
Ns is the number of nodes on the free surface S, kj and uj are the nodal values of the Lagrangian parameters and the mesh velocity at the node j. The vectors Nj 2 RNdof , where Ndof is the total number of degrees of freedom of the system. These vectors are constructed in order to add the penalization terms to the appropriate equations in the linear system. For each node j, these vectors are de®ned by Nj 0; . . . ; 0; n1j ; n2j ; n3j ; 0; . . . ; 0 where nkj are the components of the unit normal vector nj at the node j. In taking the derivation of the constraint F u with respect to the mesh velocity u, we assume that the normal vector n is constant over the time step. 4.2. Free surface discretization The nonlinear system necessitates the values of the unit normal vector at the nodes on the free surface. The nodal value of the unit normal is not well de®ned when using linear interpolation, dierent approaches for de®ning the normal vector can be used. Mass conservation in the Navier±Stokes problem is very sensitive to the particular approach used for the normal. For this reason, it constitutes one of the major problems in the ALE description. In order to determine the best approach, we consider a situation, which arises in sloshing tank problems. In these problems, ¯uid mass can only pass through the free surface boundary, and thus the mass conservation is expressed by I v n ds 0: 29 S
For an incompressible ¯uid, the conservation of mass is equivalent to the conservation of the ¯uid volume. To conserve the volume of the ¯uid domain, the mesh velocity u should satisfy
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
459
Z S
u n ds 0:
30
The free surface condition is satis®ed on each node on the free surface through the augmented Lagrangian formulation uj
vj nj 0;
j 1; . . . ; Ns :
The normal vector nj at the node j is de®ned by I nj wj n ds;
31
32
wj is the shape function at the node j. From Eq. (32), the normal vector is well de®ned. Substituting (32) into (31) and summing over the free surface nodes, we get I v u n ds 0; 33 S
PNs PNs where u j1 uj wj and v j1 vj wj . Volume conservation is deduced from equation of mass and Eq. (33). 4.3. Time integration In order to obtain a semi-discrete system governing the mesh motion, with stability properties, let us introduce a mass matrix M and a damping matrix C de®ned by ~ M mK;
~ C cK;
34
where the parameters m and c are positive and will be determined later. K~ oK d=od is the tangent stiness. Adding the mass and the damping terms to the nonlinear system (28), we get at each time step the following system: Man1 Cun1 K dn1
Ns X kn1 eF un1 Nj 0;
35
j1
where an1 is the mesh acceleration and un1 is the mesh velocity. To keep down the proliferation of notations, we dropped the index j in the penalization terms kn1 and un1 . The nonlinear system (35) is completed by the following ®nite dierence formulae: Dt2 1 2ban 2b an1 ; 2 can c an1 :
dn1 dn Dt un un1 un Dt 1
36
The parameters a and b are selected for numerical stability. To solve the systems (35) and (36), we use the HHT method for numerical dissipation Man1 Cuna K dna
Ns X kna eF una Nj 0; j1
where
37
460
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
una a un1 1
aun ;
dna a dn1 1
adn ;
kna a kn1 1
akn :
In order to solve the nonlinear system (37), a predictor±corrector algorithm is used. For the mesh motion, the so-called displacement formulation is implemented. The ®nal displacement of the previous time step is used as a predictor value for the current time step. The predictors for the displacement, velocity and acceleration are given by 0 dn ; dn1
u0n1 1 a0n1
c=bun Dt 1 c=2ban ; 1 1 un 1 an : b Dt 2b
38
A nonlinear algebraic problem is obtained for the acceleration an1 by eliminating dn1 and un1 from Eqs. (37) and (38). At each time step, the nonlinear algebraic problem for an1 is solved using Newton's method starting from the predictor solution. In order to decrease the nonlinear eects of the system, we assume the normal vector constant over the time step. The normal vector is not updated for the nonlinear Newton's method. Using the expressions of the mass matrix M and the damping matrix C, we get the following linear system to solve: i Given ain1 ; uin1 ; dn1 and kn . First, we update the Lagrange multiplier, using Uzawa's method kin1 kn
eF uin1 :
Second, we solve the following linear system for the corrector phase: ( ) Ns X m cabDt acDt2 K~ acDte Nj Nj Dai1 Mai Ri j1
n1
n1
na ;
39
where i Rina Cuina K dna
Ns X j1
kina F uina Nj :
40
Once the linear system (39) is solved, the corrector phase is given by i i1 ai1 n1 an1 Dan1 ;
i1 un1 uin1 cDt ai1 n1 ;
i1 i dn1 dn1 bDt ai1 n1 :
The parameters m and c in (35) are selected such that the linear homogeneous part of the system (36) governing the mesh motion, given by ~ K ma acu ad 0
41
is critically damped. The linear system (41) can be considered as single degree of freedom equation m
d2 x dx ca ax 0: 2 dt dt
42
The solution of the single Eq. (42) is critically damped if ma2 c2
4ma 0
43
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
and is given by x t x 0 exp
ac at : 2m
461
44
The second relationship between the parameters m and c is deduced by assuming that the solution (44) decreases by 90% at the ®rst step, which gives acDt exp 0:1: 45 2m From relations (44) and (45), we get the values of the parameters m and c m
aDt2 log 0:1
2
;
c 2 m=a
1=2
:
Remark. If the parameters a; b and c are chosen such that a 2 2=3; 1, b 1=4 2 an unconditionally stable second-order accurate scheme results, see [19].
a2 and c 1=2 3
a,
5. Numerical examples To determine the eectiveness of the ALE method, we apply the formulation to three dierent engineering problems. The ®rst is a ®nite element problem of wave propagation. The numerical simulation of this problem has been carried out, ®rst by Hughes et al. [3], where the results are compared to the experimental ones, see [12]. The second example involves a large amplitude sloshing tank. Numerical results are compared to those given by Liu et al. [14]. In the last example, a rotation of ¯uid in a cylindrical container about its central axis is considered. Analytical solution exists for the steady problem, which consists of a rigid body rotation. 5.1. The solitary wave problem The solitary wave traveling in a rectangular channel is de®ned as a single elevation above a surrounding undisturbed water level. The wave is neither followed nor preceded by another elevation of the surface. The wave travels without change of shape and essentially constant velocity. Shallow water theory is used for the analytical approach. This problem has been studied by Hughes et al. [3]. The wave is generated by prescribing the velocity time history of the ¯uid along the left-hand boundary of the domain. The physical wave parameters are selected in keeping with the numerical simulation in [3], and the experimental wave generation facility in Keck Laboratory at the California Institute of Technology. The analytical time history of displacement of the wave is given by cjt i gh d 1 tanh 4 ; j d c g D p1=2 ; 1=2 3g ; j D where g 1, p 0:86, L 9949:094 the length of the channel, and D 10 the depth of the channel. The time increment Dt h=c is selected such that cfl 1, h is the element size. The wave-generation function is designed to produce a solitary wave, one in which nonlinear dispersive effects are balanced so that the wave propagated without distortion. The computational domain consists of 160 elements mesh in the horizontal direction, and two elements through the depth. As indicated in [3], by using two elements through the depth, most of the dispersion emanating from the tail of the wave can be avoided. The numerical results are in good
462
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
agreement with those given in [3]. The variation of the wave volume during the 150 time steps is less than 2%, and the amplitude decay is less than 3%. This is due to the fact that the GLS formulation models accurately the ¯uid incompressibility. The fact is that we use several passes for both the ¯uid and the mesh stagger, instead of one, with the mass and volume conservation being accurately satis®ed. This clearly shows that the numerical method is not aected by any arti®cial damping eect. In order to get less amplitude decay and conservation of the wave volume, Newmark parameters for the time integration are considered, thus the parameters a 1, b 1=2, c 1=2 are used for both ¯uid and ALE equations. Numerical experience shows that the time derivative in the least-square terms in Eq. (16) should be considered for transient problems; when this term is not included, a greater decay of wave amplitude is obtained. For steady-state problems, this term does not need to be included in the ®nite element formulation. The successful application of the present approach to this problem indicates the possibility of employing this technique to a wide variety of water wave problems where analytical methods are dicult, and sometimes impossible to derive. Fig. 1 shows the wave shape after 50, 100, 150 time steps. 5.2. The sloshing tank To show the adaptivity of the ALE formulation for large mesh distortion, the analysis of large amplitude sloshing in a rectangular tank has been carried out. The problem has been studied by Huerta and Liu [4]. A schematic view of the tank and the boundary conditions are shown in Fig. 2. The width of the tank is L 0:8 m and the depth D 0:3 m. The tank is excited by a sinusoidal acceleration a A g sin xt;
46
where g; x and t are the gravity, the circular frequency and time, respectively. A is an arbitrary constant 1=2 governing the amplitude of the acceleration. The Reynolds number based on the reference velocity g D 5 is 5 10 . Gravity acts downwards and the amplitude of the motion is A 0:001. Using these parameters, experimental results show that the resonant frequency of the tank is x 0:88 Hz. The time increment Dt p=30x, and the sloshing period is T 2p=x. The ®nite element mesh consists of 441 21 21 elements. In order to conserve mass during the mesh motion, one uses the following strategy at each time step:
Fig. 1.
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
463
Fig. 2.
First, we run the ¯uid problem to convergence, which insures mass conservation I v n ds 0; S
47
Second, we run the mesh problem to convergence, and from free surface boundary conditions and mass conservation, we have I u n ds 0: 48 S
Eq. (48) indicates that at each time step, the mesh motion conserves the computational volume and thus mass conservation is satis®ed. Fig. 3 shows the vertical displacement of the nodes at the wall tank, as a function of time. The results are in good agreement with the results given in [4,21].
Fig. 3.
464
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
Fig. 4.
Fig. 5.
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
465
5.3. The rotating bucket Let us consider a ¯uid in a cylindrical container. The cylinder is rotating about its central axis, a schematic view is shown in Fig. 4. We assume that the container has been rotating at constant velocity long enough for the ¯uid to have attained rigid rotation. In the rotating frame, the ¯uid has no-slip boundary conditions at the wall and at the bottom of the cylinder. The ¯uid is subjected to gravity body force and the external force due to the centripetal acceleration f q x x xxr, where r is the distance to the central axis, and x is the rotational velocity. At steady-state, the free surface de®ned by the free surface condition (6) is a parabolic. Fig. 5 shows the graph of exact and numerical free surfaces. This example illustrates that the numerical integration parameters a 2=3, b 4=9 and c 9=4 lead to an accurate solution for the steady-state problem.
6. Conclusion In this paper, an ALE free surface formulation for viscous incompressible ¯ow has been presented. The GLS method is used for the ®nite element spatial discretization of the Navier±Stokes equations. A predictor±corrector scheme with numerical damping is adopted for time integration. For the mesh motion, elasticity type equations, with appropriate boundary conditions are used. Two types of constitutive models are used for the mesh motion. The ®rst one is a neo-Hookean model for the incompressible ¯ows, and draining tank problem. For the mesh motion, an appropriate discretization of the normal vector is used, in order to satisfy mass conservation in the computational domain. Numerical examples show that the ALE method is appropriate for steady and transient free surface problems.
References [1] W.F. Noh, A time dependent two space dimensional coupled Eulerian Lagrangian code, Meth. Comput. Phys. 3 (1964). [2] C.W. Hirt, A.A. Amsden, H.K. Cook, An arbitrary Lagrangian±Eulerian method for all ¯ow speeds, J. Comput. Phys. 14 (1974) 227±253. [3] T.J.R. Hughes, W.K. Liu, T.K. Zimmerman, Lagrangian±Eulerian ®nite element formulation for viscous ¯ows, Comput. Methods Appl. Mech. Engrg. 29 (1981) 329±349. [4] A. Huerta, W.K. Liu, Viscous ¯ows for large free surface motion, Comput. Methods Appl. Mech. Engrg. 69 (1988) 277±324. [5] T. Belytschko, D.F. Flanagan, J.M. Kennedy, Finite element method with user-controlled meshes for ¯uid±structure interactions, Comput. Methods Appl. Mech. Engrg. 33 (1982) 689±723. [6] O. Pironneau, Shape Optimization for Elliptic Problems, Springer, Berlin, 1983. [7] M. Souli, Y. Demay, A. Habbal, Finite element study for the draw resonance instability, Eur. J. Mech./Fluids 11 (6) (1992). [8] M. Souli, J.P. Zolesio, Shape derivative of discretized problems, Comput. Methods Appl. Mech. Engrg. 108 (1993). [9] B. Ramaswany, M. Kawahara, Arbitrary Lagrangian±Eulerian ®nite element method for unsteady convective incompressible viscous free surface ¯uid ¯ow, Int. J. Numer. Methods Fluids 7 (1987) 1053±1075. [10] R. Glowinski, O. Pironneau, Finite element methods for Navier±Stokes equations, Ann. Rev. Fluid Mech. 24 (1992) 167±204. [11] A.N. Brooks, T.J.R. Hughes, Streamline upwinded/Petrov±Galerkin methods for advection dominated ¯ows with particular emphasis on the incompressible Navier±Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199±259. [12] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new ®nite element formulation for computational ¯uid dynamics, Comput. Methods Appl. Mech. Engrg. 54 (1986) 223±234. [13] F. Shakib, Finite element analysis of the compressible Euler and Navier±Stokes equations, Ph.D. Thesis, Stanford University, 1988. [14] C. Johnson, J. Saranem, Diusion methods for the incompressible Euler and Navier±Stokes equations, Math. Comput. 47 (1986) 1±18. [15] T. Tezduyar, J. Liou, M. Behr, New strategy for FEM computations involving moving boundaries and interfaces, I, USMI Report 90/69. [16] A. Masud, A space-time ®nite element method for ¯uid±structure interactions, Ph.D. Thesis, Stanford University, 1993. [17] T.J.R. Hughes, L.P. Franca, M. Balestra, A new ®nite element formulation for computational ¯uid dynamics. V. Circumventing the Babuska±Brezzi condition, a stable Petrov±Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85±99.
466
M. Souli, J.P. Zolesio / Comput. Methods Appl. Mech. Engrg. 191 (2001) 451±466
[18] H.M. Hilbert, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms for structural dynamics, Engrg. Struct. Dyn. 5 (1977) 283±292. [19] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewoods Clis, NJ, 1987. [20] Z. Johan, T.J.R. Hughes, F. Shakib, A matrix free implicit iterative solver for compressible ¯ow problems. [21] D. Gorin, Ph.D., California Institute of Technology, 1979.