The Particle Finite Element Method: An efficient method to solve CFD problems with free-surfaces and breaking waves

c

TU Delft, The Netherlands, 2006

THE PARTICLE FINITE ELEMENT METHOD: AN

EFFICIENT METHOD TO SOLVE CFD PROBLEMS WITH

FREE-SURFACES AND BREAKING WAVES

S.R. Idelsohn1,2_{, E. O˜nate}1_{, J. Marti}2_{, M.A. Celigueta}1 _{and A. Limache}2 1 _{International Center for Numerical Methods in Engineering (CIMNE)}

Universidad Politecnica de Cataluna Campus Norte UPC, 08034 Barcelona, Spain

Web page: http://www.cimne.com

2 _{International Center for Computational Method in Engineering (CIMEC)} Universidad Nacional del Litoral and CONICET

G¨uemes 3450, Santa Fe, Argentina

Key words: Lagrangian formulation, fluid-structure, particle finite element.

Abstract. We present a general formulation for analysis of fluid-structure interaction

problems using the particle finite element method (PFEM). The key feature of the PFEM is the use of an updated Lagrangian description to model the motion of nodes (particles) in both the fluid and the structure domains. Nodes are viewed as material points (called particles) which can freely move and even separate from the main analysis domain rep-resenting, for instance, the effect of water drops. A mesh connects the nodes defining the discretized domain where the governing equations, expressed in an integral from, are solved as in the standard FEM. An incremental iterative scheme for the solution of the non linear transient coupled fluid-structure problem is described. Examples of application of the PFEM to solve a number of fluid-structure interaction problems involving large motions of the free surface and splashing of waves are presented.

1 INTRODUCTION

There is an increasing interest in the development of robust and efficient numerical methods for analysis of engineering problems involving the interaction of fluids and struc-tures accounting for large motions of the fluid free surface and the existence of fully or partially submerged bodies. Examples of this kind are common in ship hydrodynam-ics, off-shore structures, spill-ways in dams, free surface channel flows, liquid containers, stirring reactors, mould filling processes, etc.

particles is decoupled from that of the mesh nodes. Hence the relative velocity between mesh nodes and particles is used as the convective velocity in the momentum equations. Typical difficulties of FSI analysis using the FEM with both the Eulerian and ALE formulation include the treatment of the convective terms and the incompressibility con-straint in the fluid equations, the modelling and tracking of the free surface in the fluid, the transfer of information between the fluid and solid domains via the contact interfaces, the modelling of wave splashing, the possibility to deal with large rigid body motions of the structure within the fluid domain, the efficient updating of the finite element meshes for both the structure and the fluid, etc.

Most of these problems disappear if a Lagrangian description is used to formulate the governing equations of both the solid and the fluid domain. In the Lagrangian formulation the motion of the individual particles are followed and, consequently, nodes in a finite element mesh can be viewed as moving material points (hereforth called “particles”). Hence, the motion of the mesh discretizing the total domain (including both the fluid and solid parts) is followed during the transient solution.

In this paper we present an overview of a particular class of Lagrangian formulation developed by the authors to solve problems involving the interaction between fluids and solids in a unified manner. The method, called the particle finite element method (PFEM), treats the mesh nodes in the fluid and solid domains as particles which can freely move and even separate from the main fluid domain representing, for instance, the effect of water drops. A finite element mesh connects the nodes defining the discretized domain where the governing equations are solved in the standard FEM fashion. The PFEM is the natural evolution of recent work of the authors for the solution of FSI problems using Lagrangian finite element and meshless methods [Aubry et al. (2005); Idelsohn et al. (2003a; 2003b; 2004); O˜nate et al. (2003; 2004a,b)].

An obvious advantage of the Lagrangian formulation is that the convective terms dis-appear from the fluid equations. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes. Indeed for large mesh motions remeshing may be a frequent necessity along the time solution. We use an innovative mesh regeneration procedure blending elements of different shapes using an extended Delaunay tesselation [Idelsohn et al. (2003a; 2003c)]. Furthermore, this special polyhedral finite element needs special shape functions. In this paper, meshless finite element (MFEM) shape functions have been used [Idelsohn et al. (2003a)].

The layout of the paper is the following. In the next section the basic ideas of the PFEM are outlined. Next the basic equation for an incompressible flow using a Lagrangian description and the elastic solid using a hypo-elastic approximation are presented. Details of the treatment of the coupled FSI problem are given. The procedures for mesh generation and for identification of the free surface nodes are briefly outlined. Finally, the efficiency of the particle finite element method (PFEM) is shown in its application to a number of FSI problems involving large flow motions, surface waves, moving bodies. etc.

2 THE BASIS OF THE PARTICLE FINITE ELEMENT METHOD

Let us consider a domain containing both fluid and solid subdomains. The moving fluid particles interact with the solid boundaries thereby inducing the deformation of the solid which in turn affects the flow motion and, therefore, the problem is fully coupled.

In the PFEM approach presented here, both the fluid and the solid domains are mod-elled using an updated Lagrangian formulation. That is, all variables in the fluid and solid domains are assumed to be known in the current configuration at time t. The new set of variables in both domains are sought for in the next or updated configuration at time t+∆t (Figure 1). The finite element method (FEM) is used to solve the continuum equations in both domains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for both the fluid and solid problems in the standard FEM fashion. We note again that the nodes discretizing the fluid and solid domains are viewed as material particles which motion is tracked during the transient solution. This is useful to model the separation of fluid particles from the main fluid domain and to follow their subsequent motion as individual particles with a known density, an initial acceleration and velocity and subject to gravity forces.

It is important to note once more that each particle is a material point characterized by the density of the solid or fluid domain to which it belongs. The mass of a given domain is obtained by integrating the density at the different material points over the domain.

The quality of the numerical solution depends on the discretization chosen as in the standard FEM. Adaptive mesh refinement techniques can be used to improve the solution in zones where large motions of the fluid or the structure occur.

2.1 Basic steps of the PFEM

For clarity purposes we will define the collection or cloud of nodes (C) pertaining to the fluid and solid domains, the volume (V) defining the analysis domain for the fluid and the solid and the mesh (M) discretizing both domains.

A typical solution with the PFEM involves the following steps.

1. The starting point at each time step is the cloud of points in the fluid and solid domains. For instance n_{C denotes the cloud at time t = tn} (Figure 2).

(such as the free surface in fluids) may be severely distorted during the solution process including separation and re-entering of nodes. The Alpha Shape method [Edelsbrunner and Mucke (1999)] is used for the boundary definition (see Section 7).

3. Discretize the fluid and solid domains with a finite element mesh n_{M. In our work} we use an innovative mesh generation scheme based on the extended Delaunay tesselation (Section 6) [Idelsohn et al. (2003a; 2003b; 2004)].

4. Solve the coupled Lagrangian equations of motion for the fluid and the solid do-mains. Compute the relevant state variables in both domains at the next (updated) configuration for t + ∆t: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. An overview of the coupled FSI algorithm is given in the next section.

5. Move the mesh nodes to a new position n+1_{C where n + 1 denotes the time tn+ ∆t,} in terms of the time increment size. This step is typically a consequence of the solution process of step 4.

6. Go back to step 1 and repeat the solution process for the next time step.

3 EQUATIONS TO BE USED IN BOTH: FLUID AND SOLID MATERIAL 3.1 Constitutive equations for hypo-elastic solids

Let a material with a hypo- elastic constitutive equation like:

L(τij) = λsδijdll+ 2µsdij (1)

where τij = Jσij is the Kirchhoff stress tensor, dij = 12

∂Vi ∂xj + ∂Vj ∂xi 

is the rate of de-formation tensor, σij the Cauchy stress tensor, λs and µs the Lam´e parameters, J =

det(Fij) the Jacobean, being Fij = _∂X∂xi_j = _∂X∂Ui_j + δij, the deformation gradient tensor and

L(τij) = F ˙SFT, the Lie derivative , with Sij = Fij−1τijFij−T = Fij−1σijFij−TJ, the second

Piola-Kirchhoff stress tensor.

Dividing the strain and the stress in deviatory and the volumetric part

0

_V

_V

_V

_V

_V

_V

_t

_t

_V

_V

_V

_u

∆u

_{u ≡ u}

x

, u

x

, u

_t

UPDATED LAGRANGIAN FORMULATION

Initial configuration Current configuration

Next (updated) configuration

We seek for equilibrium at t + ∆t

Figure 1: Updated lagrangian description for a continuum containing a fluid and a solid domain This may be split as

L(τij ) = 2µsdij L
_τ ll 3dij  =  λs+ 2µs 3  dllδij (4)

The volumetric strain rate and the pressure will be written as εv = tr(dij) = dll and p = tr
−σij 3  =−σll 3 =− τll 3J (5)

Approaching the derivative in (4) by a finite time step F∆S  ∆t F T _{= 2µ} sd and tr  F ∆S 3∆tF T  =  λs+2µ₃s  εv (6)

Cloudn+1

C

Fluid node

Fixed boundary node Solid node

Initial “cloud” of nodes n

C

n _{C →}n _V Domainn

V

M

Figure 2: Sequence of steps to update a “cloud” of nodes from timen (t = t_n) to timen + 1 (t = t_n+ ∆t)

where: ˆ p(n)=−tr  F  Sn 3  FT  /J ; _pn+1 =−tr  F  Sn+1 3  FT  /J (8)

In the previous equations ˆ_σ
(n) and ˆ_p(n) represent the deviatory stress and the pressure at the beginning of the time step (n), but at the final time step configuration (n + 1).

Finally for the hypo-elastic material the constitutive relations may be written in one of the three following expressions:

σn+1_ij = ˆ_σij
(n)+2∆tµs J d  ij − pn+1δij (9) σn+1ij = ˆσij(n)+2∆tµ_J sdij +∆t_J  λs+2µs 3  εvδij (10) σn+1_ij = ˆ_σij(n)+2∆tµs J dij + ∆t J λsεvδij (11)

3.2 Constitutive relations for incompressible or quasi-incompressible fluids

For Newtonian fluid flows the Standard constitutive relations are:

where µL is the viscosity.

For quasi-incompressible flows, the volumetric train rate may be written as a function of the sound speed by:

Dp Dt
∆p ∆t =−C 2 ρεv ⇒ pn+1 = pn− ∆tC2ρεv = pn− ∆tκlεv (13)

where C is the sound speed and κl= C2ρ.

Then, the Cauchy stress tensor may be written in function of the volumetric strain rate by: σijn+1 = 2µLdij − pnδij + ∆tκlεvδij = 2µLdij − pnδij +  ∆tκl−2µ₃l  εvδij (14)

From (12–14) Newtonian constitutive relation for incompressible or near incompressible flow may be written in one of the three following manner:

σ_ijn+1 = 2µLdij − pn+1δij (15) σ_ijn+1 = 2µLdij − pnδij + ∆tκlεvδij (16) σ_ijn+1 = 2µLdij − pnδij +  ∆tκl− 2µ₃l  εvδij (17)

3.3 Unique constitutive equation for fluid and solid

In the following, only a unique constitutive equation will be used for both elastic solid and incompressible or nearly incompressible flow:

σijn+1 = ˆσij(n)+ 2µdij− ∆pn+1δij (18)

σijn+1 = ˆσij(n)+ 2µdij+ ∆tκεvδij (19)

σ_ijn+1 = ˆ_σij(n)_{+ 2µdij+ δεvδij} (20)

with the definitions for ˆ_σij(n)_{, µ, λ and κ given in Table 1.}

Eq.(18) will be used only in such cases in which all the domain or a part of the domain is totally incompressible, while Eq.(19) or (20) will be used in such cases in which all the domain may be considered as compressible or quasi-incompressible.

Then, depending of the material the following constant will be used: a) For the Fluid

µ = µl (21)

ˆ

σij
(n)≡ 0 (22)

ˆ

Solid Fluid Unique ˆ σ_ij(n) −pnδij σˆij(n) ∆tµs J µL µ ∆tλs J ∆tκl− 2µl 3 λ 1 J  λs+2µ₃s  κl κ

Table 1: Definition of physical parameters

κ = κl = C2ρl (24)

λ = λl= ∆tκεv (25)

κl =∞ (26)

means totally incompressible. b) For the solid part:

µ = ∆t J µs = ∆t J E 2(1 + ν) (27)

E Young module and ν Poison coefficient.

λ = λs= νE (1 + ν)(1 − 2ν) (28) κ = κs= 1 J  λs+2µs 3  (29)

3.4 Momentum conservation equations

The standard infinitesimal momentum conservation equation may be written in a La-grangian frame as:

ρai = ρDVi Dt = ∂σij ∂xj + ρfi (30) Plus εv =−∆p κ∆t (31)

And boundary conditions:

Vi = ¯Vi in ΓV (33)

It must be noted that the term DVi

Dt , represents the time material derivative (lagrangian)

of the velocity DVi

Dt = V

n+1 i −Vin

∆t where Vin+1 represents the velocity at time (n + 1) in the

position (n + 1). The convective term, normally included in fluid mechanics, are not explicitly present in this lagrangian formulation.

A weighted residual method will be used to approximate these equations:  V Wi  ρDVi Dt − ∂σij ∂xj − ρfi  dV +  Γσ Wi(σni− ¯σni)dΓ = 0 (34)  V Wp  −εv − ∆p κ∆t  dV = 0 (35) Or in weak form:  V Wiρ DVi Dt + ∂Wi ∂xj σij − WiρfidV −  Γσ Wiσ¯nidΓ = 0 (36)  V Wp  −εv− ∆p κ∆t   Γσ WidV = 0 (37) 3.5 Temporal integration Calling: an+1/2_i = V n+1 i − Vin ∆t = an+1_i + an_i 2 ⇒ a n+1 i = 2  V_in+1− V_in ∆t  − an i (38)

Then, at time (n + 1) the weak form of the weighted residual equation becomes:  V  Wi2ρ  V_in+1− V_in ∆t  − Wiρani +∂W_∂xi j σ n+1 ij − Wiρfin+1  dV − −  Γσ Wi¯σnidΓ = 0 (39)  V Wp  −εn+1 v − _κ∆t∆p  dV = 0 (40)

3.6 Case in which all the domain may be considered as compressible or quasi-incompressible

In this case, Eq.(25) is not used, and the constitutive equations used are the described in Eqs.(19) or (20).  V  Wiρ  Vi− Vin ∆t  − Wiρani + ∂W_∂xi j (ˆσ (n) ij + 2µdij + λεvδij  −  Γσ ¯ σnidΓ = 0 (43) or  V  Wi(2ρ)Vi+ ∆t∂Wi ∂xj 2µdij + ∆t ∂Wi ∂xi λεv  dV =  V [Wi(2ρ)V n i + ∆tWiρani − ∆tWiρfi] dV + ∆t  Γσ ¯ σnidΓ (44) Spatial discretization

Each of the velocity components will be interpolated using MFEM shape function [Idelsohn et al., 2003a] as:

Vi = NTQi (45)

where NT are the MFEM shape functions and Qi a vector containing the nodal value of

the velocity components.

For Galerkin residual approximations, the arbitrary weighting functions Wi are:

[W1, W2, W3] =  N0 _N0 00 0 0 _N   (46)

The equation to be solved in matrix form becomes:

3.7 Case in which all the domain or a part of it must be considered as incompressible

In this case, the only possibility is to use a mixed formulation Eq.(18) using the velocity and the pressure as unknown variables. The weighted residual equation remains:

 V  Wi2ρ  Vi− Vin ∆t  − Wiρani +∂W_∂xi j (ˆσ (n) ij + 2µdij − ∆pδij)− Wiρfi  dV − −  Γσ ¯ σnidΓ = 0 (49)  V Wp  −εv− ∆p ∆t  dV = 0 (50) or also  V  Wi(2ρ)Vi+ ∆t∂Wi ∂xj 2µd  ij − ∆t∂W_∂xi i ∆p  dV = 0  V  Wi(2ρ)Vin+ ∆tWiρani − ∆t∂W_∂xi j σˆ (n) ij + ∆tWiρfi  dV + ∆t  Γσ ¯ σnidΓ −  V  ∆tWp∂Vi ∂xi + Wp ∆p κ  dV = 0 (51)

And taking into account the definition of the deviatory strain rate:  V  Wi(2ρ)Vi+ ∆t∂Wi ∂xj µ ∂Vi ∂xi + ∆ ∂Wi ∂xj µ ∂Vj ∂xi − ∆t ∂Wi ∂xj 2µ 3 εvδij − ∆t ∂Wi ∂xi ∆p  dV =  V  Wi(2ρ)Vin+ ∆tWiρani − ∆t∂W_∂xi j σˆ (n) ij + ∆tWiρfi  dV + ∆t  Γσ ¯ σnidΓ −  V  ∆tWp∂Vi ∂xi + Wp ∆p κ  dV = 0 (52) Spatial discretization

Both, the velocity components and the pressure increment will be discretized by

Vi = NTQi and ∆p = Np∆P (53)

with Kij = Kij1 + Kij2 + Kij3 Mii=_V N(2ρ)NTdV ; Kii1 =  V ∂N ∂xj(µ) ∂NT ∂xj dV ; _Mpp=  V NpN T p dV Mij = Kij1 = 0∀i = j K2 ij =  V ∂N ∂xj(µ) ∂NT ∂xi dV Kij3 =  V ∂N ∂xi(− 2µ 3 ) ∂NT ∂xj dV Bip =_V ∂N ∂xiN T pdV ; Bpi=  V Np ∂NT ∂xi dV Gni =  V  N(2ρ)vni + N(∆tρ)an− ∆t_∂x∂N jσˆ n ij + N∆tρfi  dV +Γ_σN∆t¯σnidΓ (55)

It must be noted that this equation must be stabilized in order to avoid wiggles in the pressure solution due to the lack of the Babushka-Brezzi conditions. In this paper a Finite Calculus (FIC) formulation will be used to stabilize the solution [O˜nate (1998); 2000; 2004); O˜nate et al. (2001)].

At the end of each time step, the Cauchy stress σijn+1 are evaluated and saved for

the next time step. At the beginning of each time step, the previous Cauchy stress are considered as the second Piola-Kirchhoff stress for the present step and the are evaluated by

S_ijn ⇐ σ_ijn+1 ˆ

σn = (F SnF)/J (56)

4 TREATMENT OF CONTACT BETWEEN THE FLUID AND A FIXED BOUNDARY

The motion of the solid is governed by the action of the fluid flow forces induced by the pressure and the viscous stresses acting at the fixed boundary, as mentioned above.

The condition of prescribed velocities at the fixed boundaries in the PFEM are ap-plied in strong form to the boundary nodes. These nodes might belong to fixed external boundaries or to moving boundaries linked to the interacting solids. Contact between the fluid particles and the fixed boundaries is accounted for by the incompressibility condi-tion which naturally prevents the penetracondi-tion of the fluid nodes into the solid boundaries (Figure 3). This simple way to treat the fluid-wall contact is another attractive feature of the PFEM formulation.

5 GENERATION OF A NEW MESH

Fluid fixed boundary elementh < hcrit Fixed boundary n _C n _Γ n _V e hb> hcrit e n+1 _C n+1 _V e h > hcrit e n+1 _Γ e

Contact between fluid and fixed boundary

Fluid Air Fluid Air Air Fixed boundary

Fluid fixed boundary element

he_{< hc} n+1 _V n+1 Γ Fixed boundary u t i i t_Ve t+∆t_Ve t_Ve ₌ t+∆t_Ve

Node imoves in tangential direction due to

incompressibility This prevents the node

to penetrate into the fixed boundary

Air Fluid

There is no need for a contact search algorithm !! Contact is detected by during mesh generation !

Figure 3: Automatic treatment of contact condition at the fluid-wall interface

Figure 4: Generation of non standard meshes combining different polygons (in 2D) and polyhedra (in 3D) using the extended Delaunay technique.

of the mesh generation procedure and the derivation of the MFEM shape functions can be found in [Idelsohn et al. (2003a; 2003c; 2004)].

Once the new mesh has been generated the numerical solution is found at each time step using the fractional step algorithm described in the previous section.

6 IDENTIFICATION OF BOUNDARY SURFACES

One of the main tasks in the PFEM is the correct definition of the boundary domain. Sometimes, boundary nodes are explicitly identified differently from internal nodes. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes.

The extended Delaunay partition makes it easier to recognize boundary nodes. Con-sidering that the nodes follow a variable h(x) distribution, where h(x) is typically the minimum distance between two nodes, the following criterion has been used. All nodes on an empty sphere with a radius greater than αh, are considered as boundary nodes. In practice α is a parameter close to, but greater than one. This criterion is coincident with the Alpha Shape concept [Edelsbrunner and Mucke (1999)]. Figure 7 shows an example of the boundary recognition using the Alpha Shape technique.

Once a decision has been made concerning which nodes are on the boundaries, the boundary surface is defined by all the polyhedral surfaces (or polygons in 2D) having all their nodes on the boundary and belonging to just one polyhedron.

The correct boundary surface is important to define the normal to the surface. More-over, in weak forms (Galerkin) such as those used here a correct evaluation of the volume domain is important. In the criterion proposed above, the error in the boundary surface definition is proportional to h which is an acceptable error.

the pressure is fixed to the atmospheric value (Figure 5). We recall that each particle is a material point characterized by the density of the solid or fluid domain to which it belongs. The mass which is lost when a boundary element is eliminated due to departure of a node (a particle) from the main analysis domain is again regained when the “flying” node falls donw and a new boundary element is created by the Alpha Shape algorithm when the lost node is at a distance less than αh from the boundary.

A practical application of the method for identifying free surface particles is shown in Figure 6. The example corresponds to the analysis of the motion of a fluid within an oscillating ellipsoidal container.

Figure 5: Identification of individual particles (or a group of particles) starting from a given collection of nodes.

Figure 6: Motion of a liquid within an oscillating container. Position of the liquid particles at two different times.

6.1 Contact between the solid-solid interfaces

treated as fluid elements. Otherwise the elements are treated as contact elements where a relationship between the tangential and normal forces and the corresponding displacement is introduced so as to model elastic and frictional contact in the normal and tangential directions, respectively (Figure 7).

This algorithm has proven to be very effective and it allows to identify and model com-plex frictional contact conditions between two or more interacting solids in an extremely simple manner. Of course the accuracy on this contact model depends of the critical distance above mentioned.

Fluid domain Fixed boundary Solid t _M F_ti= -βF_viSign(V_ti) F_vi= K₁(h_c- h) – K₂V_niSign(V_ni) Fti F_ni e i Vni V_ti h < hc

Contact between solid boundaries

Contact elements are introduced

between the solid-solid interface

during mesh generation

Contact forces

Contact elements at the fixed boundary

t+∆t_M

h < hc

Solid

Solid

Contact interface

Figure 7: Contact conditions at a solid-solid interface

7 EXAMPLES

The examples chosen show the applicability of the PFEM to solve problems involving large fluid motions and FSI situations.

7.1 Rigid objects folling into water

The analysis of the motion of submerged or floating objects in water is of great interest in many areas of harbour and coastal engineering and naval architecture among others.

Figure 8: 2D simulation of the penetration and evolution of a cube and a cylinder in a water container. The colours denote the different sizes of the elements at several times.

7.2 Dragging of objects by water streams

Figure 10 shows the effect of a wave impacting on a cube representing a vehicle. This situation is typical in flooding and Tsunami situations.

7.3 Impact of sea waves on breakwaters and piers

The collapse of a water column is calculated using the present formulation and results compared with experiment results obtained by S. Koshizuka in reference [Koshizuka et all. (1995)]. Figure 11 show the experimental and the numerical results at different characteristic time step.

At time 0.1 sec, the right surfaces of the water start the disturbance due to the obstacle. At time 0.2 sec, the water surface is completely disturbed by the obstacle. The results compare considerably well with the experimental results. At 0.3 sec. collapsed water crashes to the right wall. At 0.4 sec, the water goes up along the right wall with separations and several drops. Finally, at 1,0 sec, the water along the right wall falls down and a new breaking wave will soon occur on the left wall.

Figure 12 represents the same problem but with an elastic obstacle. A unique formu-lation for incompressible flow and for the elastic obstacle has been used. No experimental results have been found for this example, but the shape of the deformation of the elastic beam as well as the free surface perturbation seems to be in agreement with the physics of the problem.

Figure 9: Detail of element sizes during the motion of a rigid cylinder within a water container.

Figures 14 and 15 evidence the potential of the PFEM to solve 3D problems involving complex interactions between water and solid objects. Figure 14 shows the simulation of the falling of two tetrapods in a water container. Finally, Figure 15 shows the motion of a collection of ten tetrapods placed in a slope under an incident wave.

Figure 14: Motion of two tetrapods falling in a water container.

8 CONCLUSIONS

Figure 15: Motion of ten tetrapods on a slope under an incident wave.

Acknowledgements

Thanks are given to Dr. F. Del Pin and Dr. N. Calvo for many useful suggestions.

REFERENCES

[1] R. Aubry, S.R. Idelsohn, E. O˜nate. Particle finite element method in fluid mechanics including thermal convection-diffusion. Computer & Structures , 83 (17-18), 1459– 1475, (2005).

[2] R. Codina, O.C. Zienkiewicz. CBS versus GLS stabilization of the incompressible Navier-Stokes equations and the role of the time step as stabilization parameter. Communications in Numerical Methods in Engineering, 18 (2), 99–112, (2002). [3] J. Donea, A. Huerta. Finite element method for flow problems. J. Wiley, (2003). [4] H. Edelsbrunner, E.P. Mucke. Three dimensional alpha shapes. ACM Trans.

[5] J. Garc´ıa, E. O˜nate. An unstructured finite element solver for ship hydrodynamic problems. J. Appl. Mech., 70, 18–26 January, (2003).

[6] S.R. Idelsohn, E. O˜nate, F. Del Pin, N. Calvo. Lagrangian formulation: the only way to solve some free-surface fluid mechanics problems. Fith World Congress on Computational Mechanics, Mang HA, Rammerstorfer FG, Eberhardsteiner J (eds), July 7–12, Viena, Austria, (2002).

[7] S.R. Idelsohn, E. O˜nate, N. Calvo, F. Del Pin. The meshless finite element method. Int. J. Num. Meth. Engng. 58 (6), 893–912, (2003a).

[8] S.R. Idelsohn, E. O˜nate, F. Del Pin. A lagrangian meshless finite element method applied to fluid-structure interaction problems. Computer and Structures, 81, 655– 671, (2003b).

[9] S.R. Idelsohn, N. Calvo, E. O˜nate. Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng. 192 (22-24), 2649–2668, (2003c).

[10] S.R. Idelsohn, E. O˜nate, F. Del Pin. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int. J. Num. Meth. Engng., 61, 964-989, (2004).

[11] S. Koshizuka, H. Tamko, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation. Comp. Fluid Dynamic Journal, 4, (1), 29–46, (1995). [12] E. O˜nate. Derivation of stabilized equations for advective-diffusive transport and

fluid flow problems. Comput. Meth. Appl. Mech. Engng., 151, 233–267, (1998). [13] E. O˜nate. A stabilized finite element method for incompressible viscous flows using a

finite increment calculus formulation. Comp. Meth. Appl. Mech. Engng., 182(1–2), 355–370, (2000).

[14] E. O˜nate. Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng., 60(1), 255–281, (2004).

[15] E. O˜nate, S.R. Idelsohn. A mesh free finite point method for advective-diffusive transport and fluid flow problems. Computational Mechanics, 21, 283–292, (1998). [16] E. O˜nate, J. Garc´ıa. A finite element method for fluid-structure interaction with

surface waves using a finite calculus formulation. Comput. Meth. Appl. Mech. Engrg.,

191, 635–660, (2001).

[18] E. O˜nate, S.R. Idelsohn, F. Del Pin. Lagrangian formulation for incompressible fluids using finite calculus and the finite element method. Numerical Methods for Scientific Computing Variational Problems and Applications, Y Kuznetsov, P Neittanmaki, O Pironneau (Eds.), CIMNE, Barcelona, (2003).

[19] E. O˜nate, J. Garc´ıa, S.R. Idelsohn. Ship hydrodynamics. Encyclopedia of Computa-tional Mechanics, E Stein, R de Borst, T.J.R. Hughes (Eds), J. Wiley, (2004a). [20] E. O˜nate, S.R. Idelsohn, F. Del Pin, R. Aubry. The particle finite element method.

An overview. Int. J. Comput. Methods, 1(2), 267-307, (2004b).

[21] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu. The finite element method for fluid dynamics. Elsevier, (2006).