A natural element updated Lagrangian approach for modelling fluid structure interactions

In this paper we present a novel methodology for the numerical simulation of fluid structure interactions in the presence of free surfaces. It is based on the use of the Natural Element Method (NEM) in an updated Lagrangian framework, together with the integration of the Navier-Stokes equations by employing a Galerkin-characteristics formulation. Tracking of the free-surface is made by employing shape constructors, in particular α- shapes. A theoretical description of the method is made and also some examples of the performance of the technique are included.


Introduction
The fact that meshless methods (Belytschko et al., 1994;Liu et al., 1995) do not suffer of mesh distortion opened a renewed interest in the last decade in Lagrangian formulations for some problems, being free-surface flows a typical example. Thus, it is possible to employ an updated Lagrangian strategy for the fluid domain, while employing a total or updated Lagrangian strategy for the solid. This approach is very convenient for some classes of problems, especially those involving drastic changes in the fluid domain geometry. Both domains are then formulated in similar frameworks and the coupling between them becomes more direct than in ALE formulations (see (Donea, 1983) or (Donea and Huerta, 2003)).
In this paper we describe mainly the fluid flow formulation proposed in the context of an updated Lagrangian strategy. We employ the α-shape-based Natural Element Method (α-NEM) (Cueto et al., 2002; to this end. At present the solid is assumed rigid, being prescribed its kinematics. This formulation posses some advantages, that include an exact interpolation along the boundary (Cueto et al., 2001), that allows for a standard, FE-like, treatment of the fluid-solid interface conditions. We firstly describe the bases of the α-NEM and then introduce the proposed numerical scheme for the integration of the Navier-Stokes equations. Finally, we include some examples that demonstrate the accuracy of the proposed scheme and also prove the potential of the technique.

Standard formulation
The NEM (Sukumar et al., 1998; is a Galerkin procedure based on the natural neighbor interpolation scheme, which in turn relies on the concepts of Voronoi diagrams and Delaunay triangulations (see Figure 1), to build Galerkin trial and test functions. These are defined as the Natural Neighbor coordinates of the point under consideration, that is, with respect to Figure 2, the value in the point x of the shape function associated with the node 1 is (Sibson, 1980;1981).

Figure 1. Delaunay triangulation and Voronoi diagram of a set of points
In addition, the NEM has other interesting properties such as linear consistency and smooth shape functions (C 1 everywhere except of the nodes). These functions are dependent on the position and density of nodes, leading to standard FE constant strain triangle shape functions, bilinear shape functions or rational quartic functions in different situations (see Figure 3 for a typical shape function). These properties permit an exact reproduction of linear displacement fields on the boundary of convex domains.

α-shape formulation
A slight modification of the way in which the Natural Neighbour interpolant is built was proposed to achieve linear interpolation also over non-convex boundaries (Cueto et al., 2001). This modification was based on the concept of α-shapes. These are a generalization of the concept of the convex hull of a cloud of points and are widely used in the field of scientific visualization and computational geometry to give a shape to a set of points. Alpha-shapes give shape to a cloud of points and are widely used in Computational Geometry despite having been developed quite recently. They were first introduced in two-dimensions by Edelsbrunner in 1983, and not generalised in three-dimensions until (Edelsbrunner and Mücke, 1994). An αshape is a generalisation of the convex hull of a cloud of points. It is a polytope that is not necessarily convex and that can be triangulated by a subset of the Delaunay triangulation, thereby maintaining the empty circumcircle criterion.
In what follows, we introduce the formal definition of a complete family of αshapes for a given set of points N, as in Edelsbrunner and Mücke (1994). Let N be a finite set of points in ℜ 3 and α a real where ∂ indicates the boundary of the ball or, more properly, the sphere or plane bounding b. That is, a ksimplex is α-exposed if an α-ball whose boundary passes through its defining points contains no other point of the set N. In this way, we can define a family of sets F k,α as the sets of α-exposed k-simplices for the given set N, fixed α and 0 ≤ k ≤ 2.
Based on these concepts, the α-shape of N, S α , is defined as the polytope whose boundary consists of the triangles in F 2,α , the edges in F 1,α and the points or vertices in F 0,α . As the α value decreases, the α-shape shrinks by the progressive development of cavities or holes. For this to occur, one or more α-balls can occupy the interior of a simplex. The α value clearly gives an intuitive measure of the maximum curvature in a region of the domain. The α-shape concept is also a generalisation of the convex hull since the α-shape for value α = 0 is identical to the initial set of points, i.e., S0 = N, and the α-shape for sufficiently high values of α is the convex hull of the given set.
An example of some α-shapes of the complete family for a given set of points distributed over the geometry of a human jaw can be seen in Figure 4. It has been demonstrated (Cueto et al., 2001) how the construction of the interpolant over an appropriate α-shape of the domain gives rise to an exact imposition of essential boundary conditions over any kind of domain (convex or not.) In addition, it enables us to track the flow front position accurately.

A natural neighbour updated Lagrangian strategy for the fluid domain
In this section we review the time integration scheme developed in (Gonzalez et al., 2006), that will be applied in the integration of the fluid flow equations. It is based on a Galerkin-characteristics formulation of the Navier-Stokes equations.

Governing equations
We consider here the problem of Fluid Dynamics at moderate Reynolds number. Thus, the governing equations can be set as follows. Consider a fluid in a region Ω of the space 2 or 3 . The fluid flow is governed by the following momentum and mass balance equations: where v represents the fluid velocity, σ the stress tensor, ρ represents fluid density and b the volumetric forces acting on the fluid. The constitutive equation for a newtonian fluid is given by: where D is the strain rate tensor, p the pressure and µ the dynamic viscosity of the fluid. To solve the problem we must prescribe an initial state as well as boundary conditions, as usual.

Time discretisation
The motion equations can be grouped to The weak form of the problem associated to Equations [5], [6] and [7] is: The second term in the right-hand side of Equation [8] represents the inertia effects. Time discretization of this term represents the discretization of the material derivative along the nodal trajectories, which are precisely the characteristic lines related to the advection operator. Thus, assuming known the flow kinematics at time t n t n ∆ − = − ) 1 ( 1 , we proceed as follows: where X represents the position at time t n-1 occupied by the particle located at position x at present time tn, i.e.: So we arrive to: where we have dropped the superscript in all the variables corresponding to the current time step.

Algorithmical issues
The most difficult term in Equation [11] is the second term of the right-hand side. The numerical integration of this term depends on the quadrature scheme employed. If we employ traditional Gauss-based quadratures on the Delaunay triangles, it will be necessary to find the position at time tn-1 of the point occupying at time tn the position of the integration point ξ k : where ω k represent the weight associated to integration point ξ k , and Ξ k corresponds to the position occupied at time tn-1 by the quadrature point ξ k . If we employ some type of nodal integration, as in (Chen et al., 2001;Gonzalez et al., 2004a), this procedure becomes straightforward, with the only need to store nodal velocities at time step t n-1 .
We discuss here the procedure to follow when employing Gauss quadratures on the Delaunay triangles. We proceed iteratively. Denoting by i the current iteration, ( 1 i ≥ ), we apply within a prescribed tolerance.
We have assumed that the number of natural neighbours of a given integration point does not change during a time step, thus needing the storage of nodal velocities at time t-1 only. It can occur that some of the nodes neighbouring the integration point at time t were not actually its neighbours at time t-1, but this does not constitute a problem, since the number of natural neighbours of a point is usually high (much bigger than three), so the quality of the interpolation is thus guaranteed. In fact, this procedure has shown to converge at a high speed, with no more than 3 iterations, at least for reasonable time steps.

Broken dam problem
The broken dam problem is classic when testing the performance of integration methods for free surface flows. We consider a rectangular column of water, initially retained by a door that is instantaneously removed at time t=0 (see Figure 5). When the door is removed, water flows under the action of gravity, considered as 9.81 m/s 2 . Density of water is 10 3 kg/m 3 , and a viscosity of 0.1 Pa·s was assumed as in other numerical simulations performed using different numerical strategies (see Duchemin et al., 2002 and references therein, for instance). The discrete model was composed of 3 364 nodes. No remeshing, addition or deletion of nodes was performed throughout the computation.  Figure 6 shows a comparison between numerical results and experimental ones, obtained from the literature (Martin and Moyce, 1952). As can be noticed, an excellent agreement was found between experimental and numerical results, despite the distortion of the triangulation. In Figure 7 the error in mass conservation is depicted, which remained always below 3%. The influence of the relationship between the parameter α and the nodal parameter h on this error was deeply analyzed in .

Water mill
In this example we study the flow generated by the movement of a water mill. The geometry of the container and the dimensions of the mill are shown in Figure 8a. The model is composed of 4 698 nodes, distributed uniformly at the initial time in a square domain of dimension 20×20cm. The sail is 10 cm long, with unit thickness. The sail rotates with constant angular speed of 0.5 rad/s. Stick boundary conditions were assumed on the reservoir walls, being the upper water surface a free boundary which evolves slightly during the simulation as noticed in Figure 8c. Time increment was set to 0.005s. being the fluid viscosity of 0.01 Pa.s.
The ability of the proposed method for describing flows in the framework of updated Lagrangian description is then fully proved.

Water mill partially submerged
The proposed method seems particularly well adapted for dealing with freesurface flows. If the sail is only partially submerged, then large-amplitude waves are expected, justifying the interest of the present simulation. For this purpose, we consider the same geometry as in the previous example, but maintaining the sail only partially submerged, as shown in Figure 9. Material parameters were chosen as in the previous example. In this case the time step was set to 0.03s. This test can be found in other references, see, for instance (Idelsohn et al., 2004). Note the appearance of a large amplitude wave on the free surface of the liquid. The geometry of the fluid and the eventual generation of drop and jets can be accurately described by the α-shapes. A deep study on this topic has been recently presented in , Gonzalez et al., 2006.

Conclusions
This paper proposes a Galerkin-characteristics updated-Lagrangian fluid flow formulation for simulation of fluid structure interaction problems. The fields approximation is based on the use of the Natural Element Method which makes it possible to work with the same cloud of nodes which moves with the material velocity, avoiding remeshing stages. The use of Lagrangian descriptions in both the solid and fluid domains greatly simplifies the formulation and numerical resolution of fluid structure interaction problems, especially those involving free-surfaces.
The application of the proposed scheme to real FSI problems (i.e., those in which the movement of the solid is coupled with the fluid one) is currently the aim of our research.