Latest revision as of 13:18, 23 May 2019

Published in Marine 2011, IV Int. Conference on Computational Methods in Marine Engineering. Selected papers, L. Eça, E. Oñate, J. García-Espinosa, T. Kvamsdal and P. Bergan (Eds.), pp. 65-81, Springer, 2013
DOI: 10.1007/978-94-007-6143-8_4

Abstract

We present some developments in the Particle Finite Element Method (PFEM) for the solution of complex coupled problems in marine, naval and harbour engineering involving fluid-soil-structure interaction (FSSI). The PFEM uses an updated Lagrangian description to model the motion of nodes (particles) in a continuum domain containing fluid, soil/rock and structures subdomains. A mesh connects the nodes defining the discretized domain where the governing equations for each of the constituent materials are solved with the FEM. The stabilization for dealing with an incompressibility material is introduced via the finite calculus (FIC) method. An incremental iterative scheme for solving the non linear transient FSSI problem is described. The procedure to model frictional contact conditions and material erosion at fluid-solid and solid-solid interfaces is described. We present examples of application of the PFEM to solve FSSI problems in marine, naval and harbour engineering such as the motion of rocks by water streams, the stability of breakwaters and constructions under sea waves, the sinking of ships and the collision of a ship with ice blocks.

1 Introduction

The analysis of problems involving the interaction of fluids, soil/rocks and structures is relevant in many areas of marine, naval and harbour engineering. Some representative examples are the study of off-shore and harbour structures under large waves, constructions hit by tsunamis, erosion of sea bed and landscape adjacent to sea shore, motion of ships under severe sea conditions, the simulation of the sinking of a ship, the collision of a ship with ice blocks, excavation and drilling problems in the sea in petroleum and gas engineering, etc.

The authors have developed in previous works a particular class of Lagrangian formulation for solving problems involving complex interactions between (free surface) fluids and solids. The so-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 mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM.

An advantage of the Lagrangian formulation is that the non-linear and non symmetric convective terms disappear from the fluid equations [1,2]. The difficulty is however transferred to the problem of adequately (and efficiently) moving the mesh nodes. The theory and applications of the PFEM are reported in [1-18].

The solution of a FSSI problem requires solving the momentum and mass balance equations for the underlying continuum that may include incompressible domains (as it is the case for a fluid and some soils). In our work we use a stabilized mixed FEM based on Finite Calculus (FIC) which allows us to solve incompressible continua with a linear approximation for the velocity and pressure variables [6,9,18,20-23].

The layout of the chapter is as follows. In the next section the key ideas of the PFEM are outlined. Next the basic equations for a general continuum using a Lagrangian description and the FIC formulation are schematically presented. Then an algorithm for the transient solution is briefly described. The treatment of the coupled FSSI problem and the method for mesh generation and for identification of the free surface nodes are outlined. The procedure for treating the frictional contact interaction between fluid, soil and structure interfaces is explained. We present several examples of application of the PFEM to solve FSSI problems in marine, naval and harbour engineering such as the motion of rocks by water streams, the stability of breakwaters and constructions hit by sea waves, the study of the sinking of ships and the collision of a ship with ice blocks.

2 The basis of the particle finite element method

Let us consider a domain containing both fluid and solid subdomains (the solid subdomain may include soil/rock materials and/or structural elements). 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 both the fluid and the solid domains are modelled using an updated Lagrangian formulation [24]. That is, all variables are assumed to be known in the current configuration at time ${\displaystyle t}$ . The new set of variables in both domains is sought for in the next or updated configuration at time ${\displaystyle t+\Delta t}$ . The finite element method (FEM) is used to solve the equations of continuum mechanics for each of the subdomains. Hence a mesh discretizing these domains must be generated in order to solve the governing equations for each subdomain in the standard FEM fashion [24-27].

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.

2.1 Basic steps of the PFEM

For clarity purposes we will define the collection or cloud of nodes (${\textstyle C}$) pertaining to the fluid and solid domains, the volume (${\textstyle V}$) defining the analysis domain for the fluid and the solid and the mesh (${\textstyle 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 ${\displaystyle {}^{n}C}$ denotes the cloud at time ${\displaystyle t=t_{n}}$ (Figure 1).
2. Identify the boundaries for both the fluid and solid domains defining the analysis domain ${\displaystyle {}^{n}V}$ in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution, including separation and re-entering of nodes. The Alpha Shape method [19] is used for the boundary definition.
3. Discretize the fluid and solid domains with a finite element mesh ${\displaystyle {}^{n}M.}$ We use an effect mesh generation scheme based on the extended Delaunay tesselation [3,4].
4. Solve the coupled Lagrangian equations of motion for the overall continuum. Compute the state variables in at the next (updated) configuration for ${\displaystyle t+\Delta t}$ : velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid.
5. Move the mesh nodes to a new position ${\displaystyle {}^{n+1}C}$ where ${\displaystyle n+1}$ denotes the time ${\displaystyle t_{n}+\Delta 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 for the next time step to obtain ${\displaystyle {}^{n+2}C}$ (Figure 1).

Figure 1:Sequence of steps to update a ``cloud´´ of nodes representing a domain containing a fluid and a solid part from time ${\displaystyle n}$ ( ${\displaystyle t=t_{n}}$ ) to time ${\displaystyle n+2}$ ( ${\displaystyle t=t_{n}+2\Delta t}$ )

3 FIC/FEM formulation for a Lagrangian continuum

3.1 Governing equations

The equations to be solved are the standard ones in Lagrangian continuum mechanics [24].

Momentum

Pressure-velocity relationship

In above equations ${\displaystyle v_{i}}$ is the velocity along the ${\displaystyle i}$th global (cartesian) axis, ${\displaystyle p}$ is the pressure (assumed to be positive in tension) ${\displaystyle \rho }$ and ${\displaystyle K}$ are the density and bulk modulus of the material, respectively, ${\displaystyle b_{i}}$ and ${\displaystyle {\sigma }_{ij}}$ are the body forces and the (Cauchy) stresses. Eqns (1) and (2) are completed with the constitutive relationships [12]:

Incompressible continuum

Compressible/quasi-incompressible continuum

where ${\displaystyle {\overset {\mbox{ˆ}}{\sigma }}_{ij}}$ are the component of the Cauchy stress tensor ${\displaystyle \left[{\overset {\mbox{ˆ}}{\sigma }}\right]}$

where ${\textstyle S}$ is the second Piola-Kirchhoff stress tensor, ${\textstyle F}$ is the deformation gradient tensor and ${\textstyle F}$ [24]. Parameters ${\displaystyle \mu }$ and ${\displaystyle \lambda }$ take the following values for a fluid or solid material:

Fluid ${\displaystyle \mu }$ : viscosity; ${\displaystyle \lambda =\Delta tK-{\frac {2\mu }{3}}}$

Solid ${\displaystyle \mu ={\frac {\Delta tG}{J}}}$ ; ${\displaystyle \lambda ={\frac {2G\nu \Delta t}{J(1-2\nu )}}}$ , where ${\displaystyle \nu }$ is the Poisson ratio, ${\textstyle G}$ is the shear modulus and ${\displaystyle \Delta t}$ the time increment.

In eqns (3) and (4), ${\displaystyle {\dot {\epsilon }}_{ij}}$ is the rate of deformation and ${\displaystyle {\delta }_{ij}}$ is the Kronecker delta. ${\displaystyle {}^{t}(\cdot )}$ denotes values at time ${\textstyle t}$.

Note that ${\displaystyle {}^{t}{\overset {\frown }{\sigma }}_{ij}=0}$ for a fluid in eqn.(4a), as the stresses depend on the rates of deformation only.

Indexes in eqns (1)-(4) range from ${\displaystyle i,j=1,n_{d}}$ , where ${\displaystyle n_{d}}$ is the number of space dimensions. These equations are completed with the standard boundary conditions of prescribed velocities and surface tractions in the mechanical problem [8-13].

3.2 Discretization of the equations

A key problem in the numerical solution of eqns (1)-(4) is the satisfaction of the mass balance condition for the fully incompressible case (i.e. ${\displaystyle K=\infty }$ in eqn (2)). A number of procedures to solve this problem exists in the finite element literature [25,27]. In our approach we use a stabilized formulation based in the so-called finite calculus (FIC) procedure [9,14,20-23]. The essence of this method is the solution of a modified mass balance equation which is written as

where ${\textstyle q}$ are weighting functions, ${\displaystyle \tau }$ is a stabilization parameter given by [9,14,21-23]

In the above, ${\textstyle h}$ is a characteristic length of each finite element and ${\displaystyle \vert v\vert }$ is the modulus of the velocity vector. In eqn (5) ${\displaystyle {\pi }_{i}}$ are auxiliary pressure projection variables chosen so as to ensure that the second term in eqn (5) can be interpreted as weighted sum of the residuals of the momentum equations and therefore it vanishes for the exact solution. The set of governing equations is completed by adding the following constraint equation [9].

where ${\displaystyle w_{i}}$ are arbitrary weighting functions.

The rest of the integral equations are obtained by applying the weighted residual technique to the governing eqns (1), (2), (3) and (5) and the boundary conditions [9,12,14,18,22,26].

We interpolate next in the finite element fashion the set of problem variables. For 3D problems these are the three velocities ${\displaystyle v_{i}}$ , the pressure ${\displaystyle p}$ and the three pressure gradient projections ${\displaystyle {\pi }_{i}.}$ In our work we use equal order linear interpolation for all variables over meshes of 3-noded triangles (in 2D) and 4-noded tetrahedra (in 3D). The resulting set of discretized equations using the standard Galerkin technique has the following form

Momentum

Pressure-velocity relationship

Pressure gradient projection

In eqns (8)-(10) ${\displaystyle {\overline {\left(\cdot \right)}}}$ denotes nodal variables, ${\textstyle {\overset {\cdot }{\overline {\left(\cdot \right)}}}=}$${\displaystyle {\frac {\partial }{\partial t}}{\overline {(\cdot )}}}$. The matrices and vectors are given in [14,18].

The solution in time of eqns (8)-(10) can be performed using any time integration scheme typical of the updated Lagrangian FEM [24]. A basic algorithm following the conceptual process described in Section 2 is presented in Box I.

Box I. Basic PFEM algorithm for a Lagrangian continuum

4 Generation of a new mesh

A key point for the success of the PFEM 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 [4].

The CPU time required for meshing grows linearly with the number of nodes. The CPU time for solving the equations exceeds that required for meshing as the number of nodes increases. As a general rule for large 3D problems meshing consumes around ${\textstyle 15\%}$ of the total CPU time per time step, while the solution of the equations (with typically 3 iterations per time step) and the system assembly consume approximately ${\textstyle 70\%}$ and ${\textstyle 15\%}$ of the CPU time per time step, respectively. These figures refer to analyses in a single processor Pentium IV PC and prove that the generation of the mesh has an acceptable cost in the PFEM. Indeed considerable speed can be gained using parallel computing techniques.

5 Identification of boundary surfaces

One of the main tasks in the PFEM is the correct definition of the boundary domain. Boundary nodes are sometimes explicitly identified. In other cases, the total set of nodes is the only information available and the algorithm must recognize the boundary nodes (Figure 2).

In our work we use an extended Delaunay partition for recognizing boundary nodes [4]. Considering that the nodes follow a variable ${\displaystyle h(x)}$ distribution, where ${\displaystyle h(x)}$ is typically the minimum distance between two nodes. All nodes on an empty sphere with a radius greater than ${\displaystyle \alpha h}$, are considered as boundary nodes. In practice ${\displaystyle \alpha }$ is a parameter close to, but greater than one. Values of ${\displaystyle \alpha }$ ranging around 1.3 have been found to be optimal in all examples analyzed. This criterion is coincident with the Alpha Shape concept [19].