Published in Comput. Methods Appl. Mech. Engrg. Vol. 197, pp. 1762–1776, 2008
doi: 10.1016/j.cma.2007.06.004
We present a general Lagrangian formulation for treating elastic solids and quasi/fully incompressible fluids in a unified form. The formulation allows to treat solid and fluid subdomains in a unified manner in fluidstructure interaction (FSI) situations. In our work the FSI problem is solved via the Particle Finite Element Method (PFEM). The PFEM is an effective technique for modeling complex interactions between floating and submerged bodies and free surface flows, accounting for splashing of waves, large motions of the bodies and frictional contact conditions. Applications of the unified Lagrangian formulation to a number of FSI problems are given.
Keywords: Lagrangian formulation, fluidstructure, particle finite element method
Fluidstructure interaction (FSI) problems are typically solved with a staggered scheme [23] by which the relevant variables at the fluid and solid subdomains are separately and sequentially (and iteratively) solved using as boundary conditions at the common fluidsolid boundaries the velocities (for the fluid problem) and the surface tractions (for the solid problems). The staggered scheme is ideal for using existing finite element codes, initially developed for fluid dynamics and solid mechanics problems, and the computing effect is mainly focused on the interfacing of the relevant data between the common fluid and solid boundaries. Indeed the FSI problem typically involves the motion of mesh nodes in both the fluid and solid domains. As a consequence, an arbitrary LagrangianEulerian (ALE) formulation [24] is used to model the governing equations for the fluid, while a standard Lagrangian formulation is used for the equations of the solid part.
In our work we propose a different route for solving FSI problems. Our goal is to solve the equations for both the fluid and solid domains using a unified lagrangian formulation. This basically means that the analysis domain, containing both fluid and solid subdomains which interact with each other, is seen as a single continuum domain with different material properties assigned to each of the interacting subdomains (i.e. the fluid and solid regions). This allows to make no distinction between fluids and solids for the numerical solution and a single computer code can be used for solving the FSI problem.
The governing equations for the fluid and solid domains (in a lagrangian frame of reference) are discretized and solved with the Particle Finite Element Method (PFEM). The PFEM treats the mesh nodes in the fluid and solid domains as moving material points (henceforth called 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 motion of the mesh discretizing the total domain (including both the fluid and solid parts) is followed during the transient solution. 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 [1,7,8,10,18,19,20]
In summary, the key ingredients of the unified formulation presented in this paper are: a) the definition of a unified constitutive equation for the fluid and solid materials, b) the use of a Lagrangian description to model the kinematic of both fluid and solid domains, and c) the use of the Particle Finite Element Method (PFEM) for redefinition of the domain boundaries and the treatment of frictional contact conditions.
An obvious advantage of the Lagrangian formulation is that the convective terms disappear 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 [7,9]. Furthermore, this special polyhedral finite element needs special shape functions. In this paper, meshless finite element (MFEM) shape functions have been used [7]. Another possibility is the use of a standard mesh generator with sliver elimination [25]. Nevertheless, standard mesh generator are som times too expensive in computational cost. An efficient mesh generator as such presented in ref. [9] is a key issue in a Lagrangian formulation.
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 hypoelastic 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.
Let us 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.
Figure 1: Sequence of steps to update a “cloud” of nodes from time () to time () 
Let a material with a hypo elastic constitutive equation like:

(1) 
where is the Kirchhoff stress tensor, is the rate of deformation tensor, the velocity along the th axis, the Cauchy stress tensor, and the Lamé parameters, the Jacobian, being the deformation gradient tensor, the th displacement component and , the time Lie derivative , with , the second PiolaKirchhoff stress tensor.
Dividing the strain and the stress in their deviatoric () and the volumetric parts

(2) 
Then

(3) 
This may be split as

(4) 
The volumetric strain rate and the pressure will be written as

(5) 
Approaching the derivative in (4) by a finite time step

(6) 
which may be written as a function of the Cauchy stress:

(7) 
where:

(8) 
In the previous equations and represent the deviatoric stress and the pressure at the beginning of the time step , but referred to the final time step configuration .
Finally for the hypoelastic material the constitutive relations may be written in terms of the following three expressions:

For Newtonian fluid flows the standard constitutive relations are:

(12) 
where is the viscosity.
For quasiincompressible flows, the volumetric train rate may be written as a strain function of the sound speed by:

(13) 
where is the sound speed and .
Then, the Cauchy stress tensor may be written in function of the volumetric strain rate by:

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

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

with the definitions for , , and given in Table 1.
Solid  Fluid  Unique 
 
 
 

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 quasiincompressible.
Then, depending of the material the following definitions will be used:
a) For the fluid

(21) 

(22) 

(23) 

(24) 

(25) 
Totally incompressible flows means and then .
b) For the solid part:

(26) 
where is the Young modulus, the Poison coefficient and and the Lamé parameters.

(27) 

(28) 
The standard infinitesimal momentum conservation equation may be written in a Lagrangian frame as:

(29) 
The equations are completed with the pressurevolumetric strain equations

(30) 
and the boundary conditions:

(31) 

(32) 
It must be noted that the term , represents the time material derivative (lagrangian) of the velocity where represents the velocity at time in the position . The convective term, normally included in the fluid mechanics equations, are not explicitly present in the lagrangian formulation.
A weighted residual method will be used to approximate above equations:

In weak form:

Let

(37) 
Then, at time the weak form of the weighted residual equation becomes:

In the following, the upper index will be dropped resulting in:

(40) 

(41) 
In this case, the constitutive equations are the described in Eqs.(19) or (20). The pressure is not a state variable:

(42) 
or

(43) 
Each of the velocity components is interpolated using the MFEM shape functions as:

(44) 
where are the MFEM shape functions and a vector containing the nodal values of the velocity components.
For Galerkin residual approximations, the arbitrary weighting functions are:

(45) 
The equation to be solved in matrix form becomes:

(46) 
with

(47) 
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:

or also

Taking into account the definition of the deviatoric strain rate:

Both the velocity components and the pressure increment are discretized by

(54) 
The equation to be solved in matrix form becomes:

(55) 
with

(56) 
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 BabuskaBrezzi conditions. In this paper a Finite Calculus (FIC) formulation will be used to stabilize the solution [12,13,14,16].
At the end of each time step, the Cauchy stress 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 PiolaKirchhoff stress for the present step and they are evaluated by

(57) 
where is the deformation gradient tensor and the defined in Section 3.1.
One of the key points for the success of the Lagrangian formulation described here is the fast regeneration of a mesh at every time step on the basis of the position of the nodes in the space domain. In our work the mesh is generated using the so called extended Delaunay tesselation (EDT) presented in [7–10]. The EDT allows one to generate non standard meshes combining elements of arbitrary polyhedrical shapes (triangles, quadrilaterals and other polygons in 2D and tetrahedra, hexahedra and arbitrary polyhedra in 3D) in a computing time of order , where is the total number of nodes in the mesh. The continuous shape functions of the elements can be simply obtained using the so called meshless finite element interpolation (MFEM). Details of the mesh generation procedure and the derivation of the MFEM shape functions can be found in [7–10].
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.
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. Considering that the nodes follow a variable distribution, where 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 , 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 [4].
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. Moreover, 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 which is an acceptable error.
The method described also allows one to identify isolated fluid particles outside the main fluid domain. These particles are treated as part of the external boundary where the pressure is fixed to the atmospheric value. 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 down and a new boundary element is created by the Alpha Shape algorithm when the lost node is at a distance less than from the boundary.
The examples chosen show the applicability of the PFEM to solve problems involving large fluid motions and FSI situations.
Only 2D quasiincompressible materials have been used in all the examples shown in the following sections.
The use of quasiincompressible formulation to approximate an incompressible flow has been largely used in the literature. The advantage of this approximation is that no stabilization is necessary to obtain smooth solutions. In our unique lagrangian formulation for both, the elastic solid and the incompressible flow, the advantage is more evident because both domains: the solid and the fluid may be solved identical in this case. The only difference between the solid and the fluid is the constitutive equation but both do not need the pressure as a state variable.
There is some type of fluidstructure interaction problems, named gravitational problems in which the introduction of a speed of sound much smaller than the real one do not disturb too much the results. The large majority of freesurface problem where the freesurface is in contact with the atmospheric pressure are inside this kind of gravitational problems and for this reason, the introduction of the real speed of sound is unnecessary and we can take advantage of this factor. The reader is referred to previous publications of the authors [10,1] to see 2D and 3D examples of FSI between totally incompressible fluid flows and rigid solids.
In all the examples shown bellow, the mesh in the solid domain is generated only once and then the nodes are moved without remeshing as in a classical finite element method. The PFEM approximation described before is only used in the fluid domain. In this way, the stresses from the previous time step remain at the element level for the next time step. Nevertheless the pressure values for the fluid domain are evaluated at the node (particles) to be transmitted to next time step with a new mesh.
The unified formulation presented has been widely used in fluid dynamics problems using PFEM [6,10,18,20]. Lecturers are asked to read previous references for PFEM validation on fluid mechanics problems. Nevertheless, for structural problems the solution with PFEM is new. For this reason a simple a pure structural dynamic example has been tested. The free vibration of a cantilever beam subjected to a suddenly applied shear stress across the other end boundary face is shown in Figure 2. A plane strain beam with length m, height , GPa, Poisson's ratio = 0.32, GPaN/cm and density =1450kg/m was discretized using 1331 equal space particles. The total shear stress was equal to 1MPa with the upper and lower faces being traction free.
Figure 2: Cantilever beam under a shear stress at the end length: Geometry 
This example was presented in [27] and compare with an analytical solution for the 1D beam theory with a correction for longitudinal shear deformation.
The maximum displacement of the beam as time function is represented in Figure 3 and compare with the analytical solution. Different time step were used in order to test de numerical diffusion of the method. Courant number equal 0.5 and 50 were tested. We can conclude that the approximation used in the solid part of the domain has a small numerical diffusion and that time step order of Courant number between 0.5 and 50 are enough to have excellent results.
Figure 3: Cantilever beam under a shear stress at the end length: Maximum displacement for Courant numbers 0.5 and 50 
The collapse of a water column illustrated in Figure 4 is calculated using the present formulation and compared with experiment results obtained by S. Koshizuka [11]. Figure 5 shows 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 very well with the experimental results. At 0.3 sec. the 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 4: Initial geometry of the water column 
Figure 5: Comparison between experimental and numerical results of the collapse of a water column with an obstacle 
After 0.4 sec a large deviation between the experimental results and the numerical one is observed. The reason is that the experimental results were performed in recipient with water in contact with the external air. After 0.4 sec, an air bubble is formed with the jet of water and the recipient. This incompressible air bubble changes significantly the experimental results.
In order to improve the numerical results, the same problem was solved including the recipient water and air. Figure 6 shows the new results in which the air is represented by a fluid with their corresponding density and viscosity. The results are now in a better agreement with the experimental ones and show the powerful of the PFEM to handle large differences in fluid physical characteristic.
Figure 6: Comparison between experimental and numerical results of the collapse of a water column with an obstacle. Case with water and air. 
Figure 7 represents the same problem but with a elastic obstacle with density , Young's modulus and Poisson's ratio . The geometry of the more slender obstacle is of width and height . This example was also solved in Ref. 26 with a staggered and levelset method for the free surface flow. A unified formulation for incompressible flow and for the flexible obstacle has been used in this paper. 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.
The left upper corner of the solid first gets a defection to the left when the water acts on its lower part and moves to the right while the water rises. It obtains its maximum deflection (mark (a) in Figure 8) when the water jet passes the top and is fully attached to the left side of the solid. Finally, the impact of the fluid causes the thin solid to start oscillating (b). This oscillation is damped (c) by the water located left and right of the wall. The results obtained in Ref. 26 are also presented in Figure 8.
Figure 7: Collapse of a water column with an elastic obstacle 
Figure 8: Collapse of a water column with an elastic obstacle: History of the displacement and comparison with Ref. 26 
In order to increase the interaction between the fluid and the elastic solid, a similar problem to the previous one has been solved but using a taller elastic beam builtin on the bottom of the recipient. The obstacle with density , Young's modulus and Poisson's ratio . Geometry is illustrated in Figure 9.
Figure 9: Initial geometry of a high elastic bar under a lateral wave 
Now, the system is simulated for with a constant time step . The hypoelastic solid is deflected by the impulse of the fluid wave (Figure 10).
Figure 10: A high elastic beam builtin on the bottom of the recipient 
The analysis of the motion of submerged or floating object in water is of great interest in many areas of harbour and costal engineering and naval architecture among others. In the following figures a series of examples will be described in order to show the potential of the PFEM to handle fully coupled elastic solid and viscous fluid flows in problems in which the position of the body depends totally on the internal forces produces by the fluid. The presence of free surfaces introduces another complexity to the problem which is solved easily by PFEM.
Figure 11 shows the penetration and motion of an elastic cube in a container with water. The difference in the density of the elastic body with the water results that the body will float on the freesurface. The density for the body used was 5% bigger than fluid. The colours denote the pressure distribution at different time instants. In this case, the speed of sound used for the water was .The small value for the sound speed may be observed in some time fluctuation of the pressure values, which for this kind of gravitational problems, do not affect considerably the geometric results.
Figure 12 shows the evolution of an elastic floating plate on a freesurface with breaking waves. Here again a unified fluidsolid formulation has been used. The floating plate has a density 5% bigger than the fluid. For this reason the plate becomes totally submerged in the water.
Figure 11: Elastic cube falling into a container with water 
Figure 12: Dragging an elastic plate in a water container 
A unified Lagrangian finite element framework in conjunction with the Particle Finite Element Method is ideal to treat problems involving fluids with free surfaces and submerged or floating structures. Problems such as fluidstructure interaction, large motion of fluid or solid particles, surface waves, water splashing, separation of water drops, etc. can be easily solved with the PFEM. The success of the method lies in the accurate and efficient solution of the equations of an incompressible fluid and of solid dynamics using an updated Lagrangian formulation. Other essential solution ingredients are the efficient regeneration of the finite element mesh using an extended Delaunay tesselation, the meshless finite element interpolation (MFEM), the identification of the boundary nodes using an Alpha Shape type technique and the simple algorithm to treat contact conditions at fluidsolid and solidsolid interfaces via mesh generation. The examples presented have shown the great potential of the PFEM for solving a wide class of practical FSI problems.
The authors thank Prof. J. Oliver and Dr. R. Rossi for many useful discussions. This work was partially supported by project SEDUREC of the Consolider Programme of the Ministerio de Educación y Ciencia of Spain.
[1] R. Aubry, S.R. Idelsohn, E. Oñate. Particle finite element method in fluid mechanics including thermal convectiondiffusion. Computer & Structures , 83 (1718), 1459–1475, (2005).
[2] R. Codina, O.C. Zienkiewicz. CBS versus GLS stabilization of the incompressible NavierStokes 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. Graphics, 13, 43–72, (1999).
[5] J. García, E. Oñate. An unstructured finite element solver for ship hydrodynamic problems. J. Appl. Mech., 70, 18–26 January, (2003).
[6] S.R. Idelsohn, E. Oñate, F. Del Pin, N. Calvo. Lagrangian formulation: the only way to solve some freesurface 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ñate, 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ñate, F. Del Pin. A lagrangian meshless finite element method applied to fluidstructure interaction problems. Computer and Structures, 81, 655–671, (2003b).
[9] S.R. Idelsohn, N. Calvo, E. Oñate. Polyhedrization of an arbitrary point set. Comput. Method Appl. Mech. Engng. 192 (2224), 2649–2668, (2003c).
[10] S.R. Idelsohn, E. Oñate, F. Del Pin. The particle finite element method: a powerful tool to solve incompressible flows with freesurfaces and breaking waves. Int. J. Num. Meth. Engng., 61, 964989, (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ñate. Derivation of stabilized equations for advectivediffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Engng., 151, 233–267, (1998).
[13] E. Oñate. 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ñate. Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng., 60(1), 255–281, (2004).
[15] E. Oñate, S.R. Idelsohn. A mesh free finite point method for advectivediffusive transport and fluid flow problems. Computational Mechanics, 21, 283–292, (1998).
[16] E. Oñate, J. García. A finite element method for fluidstructure interaction with surface waves using a finite calculus formulation. Comput. Meth. Appl. Mech. Engrg., 191, 635–660, (2001).
[17] E. Oñate, C. Sacco, S.R. Idelsohn. A finite point method for incompressible flow problems. Comput. Visual. in Science, 2, 67–75, (2000).
[18] E. Oñate, 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ñate, J. García, S.R. Idelsohn. Ship hydrodynamics. Encyclopedia of Computational Mechanics, E Stein, R de Borst, T.J.R. Hughes (Eds), J. Wiley, (2004a).
[20] E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry. The particle finite element method. An overview. Int. J. Comput. Methods, 1(2), 267307, (2004b).
[21] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu. The finite element method for fluid dynamics. Elsevier, (2006).
[22] O.C. Zienkiewicz, R.L. Taylor, The finite element method for solid and structural mechanics. Elsevier, (2005).
[23] Ch. Farhat, G. Kristoffer, V. Zee, Ph. Geozine. Probably secondorder time accurate looselycoupled solution algorithm for transient nonlinear computational aeroelasticity. Comput. Meth. Appl. Mech. Engrg., 195(17–18), 1973–2001, (2006).
[24] Ch. Foster, W.A. Wall, E. Ramm. Artificial added mass instabiliting in sequential staggered coupling of nonlinear structures and incompressible flows. Comput. Meth. Appl. Mech. Engrg., 196, 1278–1293, (2007).
[25] R. Lohner. A parallel advancing front grid generation scheme. AIAA, 001005, (2000).
[26] E. Walhorn, A. Kolke, B. Hubner, D. Dinkler. Fluidstructure coupling with monolithic model involving free surface flows. Computers & Structures, 83, 2100–2111, (2005).
[27] C.J. Greenshields, H.G. Weller. A unified formulation for continuum mechanics applied to fluidstructure interaction in flexible tubes. Int. J. Num. Meth. Engng., 64, 1575–1593. (2005).
Published on 01/01/2008
DOI: 10.1016/j.cma.2007.06.004
Licence: CC BYNCSA license
Are you one of the authors of this document?