# A natural element updated Lagrangian approach for modelling fluid structure interactions

Andrés Galavis, David González, Elías Cueto, Francisco Chinesta, Manel Doblaré

<figure><img src="image_1.png" /></figure>

## To cite this version:

Andrés Galavis, David González, Elías Cueto, Francisco Chinesta, Manel Doblaré. A natural element updated Lagrangian approach for modelling fluid structure interactions. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2007, 16 (3-4), pp.323 - 336. 10.3166/remn.16.323-336. [hal-01004905](https://hal.science/hal-01004905)

HAL Id: hal-01004905

https://hal.science/hal-01004905

Submitted on 18 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A natural element updated Lagrangian
approach for modelling fluid structure
interactions

Andrés Galavís* — David González* — Elias Cueto*
Francisco Chinesta**— Manuel Doblaré*

* Aragón Institute of Engineering Research, University of Zaragoza
Edificio Betancourt, María de Luna, 7, E-50018 Zaragoza, Spain
ecueto@unizar.es

** Laboratoire de mécanique des systèmes et des procédés
UMR CNRS-ENSAM-ESEM
Ecole Nationale Supérieure des Arts et Métiers, 151 boulevard de l'hôpital
F-75013 Paris

ABSTRACT. 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.

RÉSUMÉ. Dans cet article on présente une nouvelle méthodologie pour la simulation numérique des interactions fluide-structure en présence de surfaces libres. Elle se base sur l'emploi de la méthode des éléments naturels (MEN) dans un cadre Lagrangien actualisé, utilisant en même temps l'intégration des équations de Navier-Stokes à travers la méthode des caractéristiques-Galerkin. Des constructeurs de forme sont employés pour le suivi de la surface libre, particulièrement des formes-α. On présente une description théorique de la méthode, ainsi que quelques exemples du fonctionnement de la technique.

KEYWORDS: fluid structure interaction, meshless methods, natural element method, α-shapes.

MOTS-CLÉS : fluide structure interaction, méthode sans maillage, méthode des éléments naturels, formes-α.
# 1. 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; 2003) 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.

# 2. The natural element method

## 2.1. Standard formulation

The NEM (Sukumar et al., 1998; Cueto et al., 2003) 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).

$$ \phi_1(x) = \frac{A_{abfe}}{A_{abcd}} \quad [1] $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Delaunay triangulation and Voronoi diagram of a set of points</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Definition of natural neighbour coordinates</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Typical Sibson shape function (courtesy N. Sukumar)</figcaption></figure>

**2.2. α-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 $\mathbb{R}^3$ and $\alpha$ a real number with $0 \le \alpha \le \infty$. A $k$-simplex $\sigma_T$ with $0 \le k \le 3$, is defined as the convex hull of a subset $T \subseteq N$ of size $|T| = k+1$. Let $b$ be an α-ball, that is, an open ball of radius $\alpha$. A $k$-simplex $\sigma_T$ for $0 \le k \le 2$ is said to be α-exposed if there exists an empty α-ball $b$ with $T = \partial b \cap N$, where $\partial$ indicates the boundary of the ball or, more properly, the sphere or plane bounding $b$. That is, a $k$-simplex 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,\alpha}$ as the sets of α-exposed $k$-simplices for the given set $N$, fixed $\alpha$ and $0 \le k \le 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., $S_0 = 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.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Five elements of the complete family of α-shapes of a cloud of points distributed over the geometry of a human jaw. Increasing values of α from 0 (top line, left) to infinity (bottom line, right)</figcaption></figure>

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.

### 3. 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.
### 3.1. 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 $\Omega$ of the space $\mathbb{R}^2$ or $\mathbb{R}^3$. The fluid flow is governed by the following momentum and mass balance equations:

$$ \rho(\mathbf{v}_t + (\mathbf{v} \cdot \nabla)\mathbf{v}) = \nabla \boldsymbol{\sigma} + \rho \mathbf{b} \quad \text{in } \Omega \times [0,t] \quad [2] $$

$$ \nabla \cdot \mathbf{v} = 0 \quad \text{in } \Omega \times [0,t] \quad [3] $$

where $\mathbf{v}$ represents the fluid velocity, $\boldsymbol{\sigma}$ the stress tensor, $\rho$ represents fluid density and $\mathbf{b}$ the volumetric forces acting on the fluid. The constitutive equation for a newtonian fluid is given by:

$$ \boldsymbol{\sigma} = -p\mathbf{I} + \boldsymbol{\tau} = -p\mathbf{I} + 2\mu\mathbf{D} \quad [4] $$

where $\mathbf{D}$ is the strain rate tensor, $p$ the pressure and $\mu$ the dynamic viscosity of the fluid. To solve the problem we must prescribe an initial state as well as boundary conditions, as usual.

### 3.2. Time discretisation

The motion equations can be grouped to

$$ \nabla \cdot \boldsymbol{\sigma} + \rho \mathbf{b} = \rho \frac{d\mathbf{v}}{dt} = \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} \right) \quad [5] $$

$$ \nabla \cdot \mathbf{v} = 0 \quad [6] $$

$$ \boldsymbol{\sigma} = -p\mathbf{I} + 2\mu\mathbf{D} \quad [7] $$

The weak form of the problem associated to Equations [5], [6] and [7] is:

$$ \begin{gathered} \int_{\Omega} 2\mu \mathbf{D} : \mathbf{D}^* d\Omega - \int_{\Omega} p \mathbf{I} : \mathbf{D}^* d\Omega = - \int_{\Omega} \rho \mathbf{b} \mathbf{v}^* d\Omega + \int_{\Omega} \rho \frac{d\mathbf{v}}{dt} \mathbf{v}^* d\Omega \quad [8] \\ \int_{\Omega} \nabla \cdot \mathbf{v} p^* d\Omega = 0 \end{gathered} $$

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-1} = (n-1)\Delta t$, we proceed as follows:

$$ \int_{\Omega} \rho \frac{d\mathbf{v}}{dt} \mathbf{v}^* d\Omega = \int_{\Omega} \rho \frac{\mathbf{v}^n(\mathbf{x}) - \mathbf{v}^{n-1}(\mathbf{X})}{\Delta t} \mathbf{v}^* d\Omega \quad [9] $$

where $\mathbf{X}$ represents the position at time $t_{n-1}$ occupied by the particle located at position $\mathbf{x}$ at present time $t_n$, i.e.:

$$ \mathbf{x} = \mathbf{X} + \mathbf{v}^{n-1}(\mathbf{X})\Delta t \quad [10] $$

So we arrive to:

$$ \int_{\Omega} 2\mu \mathbf{D} : \mathbf{D}^* d\Omega - \int_{\Omega} p \mathbf{I} : \mathbf{D}^* d\Omega - \int_{\Omega} \rho \frac{\mathbf{v}\mathbf{v}^*}{\Delta t} d\Omega = - \int_{\Omega} \rho \mathbf{b}\mathbf{v}^* d\Omega - \int_{\Omega} \rho \frac{\mathbf{v}^{n-1}\mathbf{v}^*}{\Delta t} d\Omega \quad [11] $$

where we have dropped the superscript in all the variables corresponding to the current time step.

### 3.3. 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 $t_{n-1}$ of the point occupying at time $t_n$ the position of the integration point $\xi_k$:

$$ \int_{\Omega} \rho \frac{\mathbf{v}^{n-1} \mathbf{v}^*}{\Delta t} d\Omega = \sum_k \rho \frac{\mathbf{v}^{n-1}(\Xi_k) - \mathbf{v}^*(\xi_k)}{\Delta t} \omega_k \quad [12] $$

where $\omega_k$ represent the weight associated to integration point $\xi_k$, and $\Xi_k$ corresponds to the position occupied at time $t_{n-1}$ by the quadrature point $\xi_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, $(i \ge 1)$, we apply
$$x_k = X_k^i + v_{n-1}(X_k^{i-1})\Delta t \text{ with } x_k = X_k^0$$

[13]

until $X_k^i \approx X_k^{i-1}$ 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.

# 4. Numerical examples

## 4.1. 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 \text{ m/s}^2$. Density of water is $10^3 \text{ kg/m}^3$, and a viscosity of $0.1 \text{ Pa} \cdot \text{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><img src="image_7.png" /><figcaption>**Figure 5.** Geometry of the broken dam problem</figcaption></figure>
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 $\alpha$ and the nodal parameter $h$ on this error was deeply analyzed in (Martinez et al., 2004).

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Front position (in non-dimensional form) in time. Numerical results vs. experimental ones</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** *Error in mass conservation*</figcaption></figure>
## 4.2. 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.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) to (c) vector plot of the velocity field at three different time steps. (d) Contour plot of the velocity field</figcaption></figure>
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.

### 4.3. Water mill partially submerged

The proposed method seems particularly well adapted for dealing with free-surface 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.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** *Initial geometry of the mill problem for the partially submerged configuration*</figcaption></figure>

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 (Martinez *et al.*, 2004, Gonzalez *et al.*, 2006).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Six snapshots of the generation of a wave during the rotation of the mill</figcaption></figure>

**5. 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.

Acknowledgement

The support from the Spanish Ministry of Education and Science, through project
DPI-2005-08727-C02-01 is gratefully acknowledged.

**6. References**

Belytschko T., Lu Y. Y., Gu L., “Element-Free Galerkin Methods”, *International Journal for Numerical Methods in Engineering*, 37, 1994, p. 229-256.

Chen Jiun-Shyan, Wu Cheng-Tang, Yoon Sangpil, You Yang, “A stabilized conforming nodal integration for galerkin mesh-free methods”, *International Journal for Numerical Methods in Engineering*, 50, 2001, p. 435-466.

Cueto E., Calvo B., Doblare M., “Modeling three-dimensional piece-wise homogeneous domains using the α-shape based Natural Element Method”, *International Journal for Numerical Methods in Engineering*, 54, 2002, p. 871-897.

Cueto E., Cegoñino J., Calvo B., Doblare M., “On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods”, *Communications in Numerical Methods in Engineering*, vol. 19, n° 5, 2003, p. 361-376.

Cueto E., Doblare M., Gracia L., “Imposing essential boundary conditions in the Natural Element Method by means of density-scaled alpha-shapes”, *International Journal for Numerical Methods in Engineering*, vol. 49, n° 4, 2000, p. 519-546.

Cueto E., Sukumar N., Calvo B., Martinez M. A., Cegoñino J., Doblare M., “Overview and recent advances in Natural Neighbour Galerkin methods”, *Archives of Computational Methods in Engineering*, vol. 10, n° 4, 2003, p. 307-384.

Donea J., “Arbitrary Lagrangian-Eulerian finite element methods”, Belytschko T., Hughes T. J. R., editors, *Computer methods for transient analyses*, Elsevier-Amsterdam, 1983, p. 474-516.

Donea J., Huerta A., *Finite Element Methods for Flow Problems*, John Wiley and sons, 2003.

Duchemin L., Popinet S., Josserand C., Zaleski S., “Jet formation in bubbles bursting at a free surface”, *Physics of Fluids*, 14, 2002, p. 3000-3008.

Edelsbrunner H., Muecke E., “Three dimensional alpha shapes”, *ACM Transactions on Graphics*, 13, p. 43-72, 1994.

Gonzalez D., Cueto E., Martinez M. A., Doblare M., “Numerical integration in Natural Neighbour Galerkin methods”, *International Journal for Numerical Methods in Engineering*, vol. 60, n° 12, 2004, p. 20771-2104.
Gonzalez D., Cueto E., Doblare M., "Volumetric locking in Natural Neighbour Galerkin methods", *International Journal for Numerical Methods in Engineering*, vol. 61, n° 4, 2004, p. 11-632.

González D., Cueto E., Chinesca F., Doblaré M., "A Natural Element updated Lagrangian strategy for free-surface Fluid Dynamics", Submitted, 2006.

Idelsohn S. R., Oñate E., del Pin F., "The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves", *International Journal for Numerical Methods in Engineering*, vol. 61, n° 7, 2004, p. 964-989.

Liu W. K., Jun S., Li S., Adee J., Belytschko T., "Reproducing kernel particle methods", *International Journal for Numerical Methods in Engineering*, 38, 1995, p. 1655-1679.

Martin J.C., Moyce W.J., "Part IV. an experimental study of the collapse of liquid columns on a rigid horizontal plane", *Phil. Tran. R. Soc. London*, 244, p. 312, 1952.

Martinez M. A., Cueto E., Alfaro I., Doblare M., Chinesta F., "Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods", *International Journal for Numerical Methods in Engineering*, vol. 60, n° 13, 2004, p. 2105-2129.

Sibson R., "A Vector Identity for the Dirichlet Tesselation", *Mathematical Proceedings of the Cambridge Philosophical Society*, 87, 1980, p. 151-155.

Sibson R., "A brief description of natural neighbour interpolation", *Interpreting Multivariate Data*, V. Barnett (Editor), John Wiley, 1981, p. 21-36.

Sukumar N., Moran B., Belytschko T., "The Natural Element Method in Solid Mechanics", *International Journal for Numerical Methods in Engineering*, vol. 43, n° 5, 1998, p. 839-887.

Sukumar N., Moran B., Yu Semenov A., Belikov V. V., "Natural Neighbor Galerkin Methods", *International Journal for Numerical Methods in Engineering*, vol. 50, n° 1, p. 1-27, 2001.