|
|
| Line 1: |
Line 1: |
| − | <!-- metadata commented in wiki content
| + | Published in ''Comput. Meth. Appl. Mech. Engng.'', Vol. 195 (17-18), pp. 2100-2123<br /> |
| | + | doi: 10.1016/j.cma.2005.02.026 |
| | | | |
| − |
| |
| − | <div id="_GoBack" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <big>'''Fluid-Structure Interaction Using the'''</big></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <big>'''Particle Finite Element Method'''</big></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <big>'''S.R. Idelsohn<sup>(1,2)</sup>, E. Oñate<sup>(2)</sup>, F. Del Pin<sup>(1)</sup>''' '''and Nestor Calvo<sup>(1)</sup>'''</big></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">''(1) ''International Center for Computational Methods in Engineering (CIMEC)</span></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">Universidad Nacional del Litoral and CONICET, Santa Fe, Argentina</span></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">e–mail: [mailto:sergio@ceride.gov.ar sergio@ceride.gov.ar]</span></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">''(2) ''International Center for Numerical Methods in Engineering (CIMNE)</span></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">Universidad Politécnica de Cataluña, Barcelona, Spain</span></div>
| |
| − |
| |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| |
| − | <span style="text-align: center; font-size: 75%;">e–mail: [mailto:onate@cimne.upc.es onate@cimne.upc.es]</span></div>
| |
| − | -->
| |
| | | | |
| | ==Abstract== | | ==Abstract== |
| | | | |
| − | ''In the present work a new approach to solve fluid-structure interaction problems is described. Both, the equations of motion for fluids and for solids have been approximated using a material (lagrangian) formulation. To approximate the partial differential equations representing the fluid motion, the shape functions introduced by the Meshless Finite Element Method (MFEM) have been used. Thus, the continuum is discretized into particles that move under body forces (gravity) and surface forces (due to the interaction with neighboring particles). All the physical properties such as density, viscosity, conductivity, etc., as well as the variables that define the temporal state such as velocity and position and also other variables like temperature are assigned to the particles and are transported with the particle motion. The so called Particle Finite Element Method (PFEM) provides a very advantageous and efficient way for solving contact and free-surface problems, highly simplifying the treatment of fluid-structure interactions. ''
| + | In the present work a new approach to solve fluid-structure interaction problems is described. Both, the equations of motion for fluids and for solids have been approximated using a material (lagrangian) formulation. To approximate the partial differential equations representing the fluid motion, the shape functions introduced by the Meshless Finite Element Method (MFEM) have been used. Thus, the continuum is discretized into particles that move under body forces (gravity) and surface forces (due to the interaction with neighboring particles). All the physical properties such as density, viscosity, conductivity, etc., as well as the variables that define the temporal state such as velocity and position and also other variables like temperature are assigned to the particles and are transported with the particle motion. The so called Particle Finite Element Method (PFEM) provides a very advantageous and efficient way for solving contact and free-surface problems, highly simplifying the treatment of fluid-structure interactions. |
| − | | + | |
| − | '''Key words: '''Fluid-Structure interaction, Particle methods, Lagrange formulations, Incompressible Fluid Flows, Meshless Methods, Finite Element Method.
| + | |
| − | | + | |
| − | == Introduction==
| + | |
| − | | + | |
| − | Many classifications have been proposed to enclose the numerical formulations that approximate the continuum equations that govern incompressible fluid flows. In particular the one describing the way that convection is treated divides the numerical formulations into two classes, namely ''material (or lagrangian) formulations'' and ''spatial (or eulerian) formulations''. The first one describes convection by placing a set of axes over the material particles that move accordingly to the equations of motion. In the eulerian case the axes are set fixed in space and convection terms are included in the equations describing the transport of the fluid flow. The present work will describe a method that uses a material formulation. The equations of motion for both, the solid and fluid do not present convection terms, implying that the convection effect is directly obtained by moving the discrete domain.
| + | |
| − | | + | |
| − | Many authors have taken advantage of lagrangian formulations to describe different types of problems. The ''Smooth Particle Hydrodynamics'' (SPH) method developed by Monaghan(1977)[mon81, mon97] should be mentioned as a pioneer method of this kind.
| + | |
| − | | + | |
| − | Many other methods have been derived from SPH. One that has shown remarkable results is the ''Moving Particle Semi-Implicit'' method (MPS) introduced by Koshizuka and Oka(1996)[kos96]. These methods use a kernel function to interpolate the unknowns. SPH uses a weak formulation while MPS uses a strong form of the governing equations.
| + | |
| − | | + | |
| − | Ramaswamy (1986)[ram87] proposed a lagrangian finite element formulation for a 2-D incompressible fluid flow. In that paper the mesh was convected according to the equations of motion but without change of topology, making it rather limiting when the elements got highly distorted. The equations of motion were discretized in space by using the finite element method with linear shape functions.
| + | |
| − | | + | |
| − | Another possible classification for numerical formulations may be the one that separates the methods that make use of a ''standard finite element mesh'' (like those made of tetrahedra or hexahedra), and the methods that do not need a standard mesh, namely the ''meshless methods''. The formulation described in this paper can be considered a particular class of meshless method. Again, SPH might be cited as one the first meshless methods.
| + | |
| − | | + | |
| − | Indeed, after Monaghans work and in particular in the past 20 years, many have been the attempts to develop a robust meshless method that could approximate PDE’s in 2-D and 3-D with acceptable accuracy, convergence and speed. Among others, the methods based on ''Moving Least Square'' interpolations [nay92,'' ''bel94], ''Partition of Unity'' [dua95], and the ones based on the ''natural neighbor'' interpolation functions [bra95, suk98] may be listed.
| + | |
| − | | + | |
| − | In this work the interpolation function used by the Meshless Finite Element Method (MFEM) [ide03a] will be implemented. This function uses the Voronoï diagram of the cloud of points to construct the interpolant. The ''extended Delaunay tessellation'' (EDT) [ide03b] is applied to connect the neighboring particles. The EDT provides polyhedral elements that are sliver-free in 3-D, avoiding instabilities of the Delaunay tessellation due to distorted tetrahedra. The MFEM shape functions adapt automatically to the polyhedra and in the case that the polyhedron is a simplex, the shape function behaves exactly as the linear finite element shape function.
| + | |
| − | | + | |
| − | Fluid-structure interaction (FSI) problems have been of special interest for designers and engineers in the past 20 years. This explains why more robust and stable formulations have been developed to assist the approximation of contact problems. Embedded methods have been developed by Löhner et al. [loh03] where a single mesh is used to partition the fluid as well as the structure. Also Arbitrary Lagrangian-Eulerian (ALE) formulations [Sou00] have given acceptable results when the displacements or the geometry deformations are not excessively large.
| + | |
| − | | + | |
| − | The approximation for the FSI problem depends basically on the coupling of the fluid and structure equations. Based on this coupling FSI problems may be divided into problems with weak interaction and problems with strong interaction. The later are found when elastic deformation of the solid takes place. The weak interpolation case happens when large rigid displacements are present. This situation is typical in ship hydrodynamics, when a rigid body moves according to the forces given by the pressure field obtained from the fluid dynamic problem. These forces applied to the rigid body will accelerate it, changing its velocity and therefore, its position.
| + | |
| − | | + | |
| − | FSI problems have been classically solved in a partitioned manner solving iteratively the discretized equations for the flow and the solid domain separately. The solution of both, fluid flow and solid, with the same material formulation, open the door to solve the global coupled problem in a monolithic fashion. Nevertheless, in this paper the rigid solid will still be solved separately from the fluid. A partitioned method [pip95, mok01] or iterative method [rug00, rug01, zha01 is chosen to solve the coupling between the fluid and solid. The advantage to use a material formulation for both, solid and fluid parts will be used here only to better reproduce breaking waves or separated drops in the fluid, which are phenomena impossible to reproduce using a spatial formulation.
| + | |
| − | | + | |
| − | The layout of the paper is the following: in the next section the basic lagrangian equations of motion for the fluid and solid domains are given. Next the discretization method chosen to solve the incompressible fluid flow equations and the solid dynamics in time equations are detailed. The algorithm for the recognition of the boundary nodes and the treatment of the free-surface in the fluid is explained. Finally the efficiency of the Particle Finite Element Method for solving a variety of fluid-structure interaction problems involving large motion of the free-surface in the fluid is shown.
| + | |
| − | | + | |
| − | ==Equations of motion==
| + | |
| − | | + | |
| − | ===2.1 Fluid dynamic problem: updating the fluid particle positions===
| + | |
| − | | + | |
| − | The fluid particle positions will be updated via solving the lagrangian form of the Navier-Stokes equations.
| + | |
| − | | + | |
| − | Let [[Image:Draft_Samper_772971023-image1.png|18px]] the initial position of a particle a time ''t=t<sub>0</sub> ''and let [[Image:Draft_Samper_772971023-image2.png|12px]] the final position. Been [[Image:Draft_Samper_772971023-image3.png|84px]] the velocity of the particle in the final position the following approximate relation can be written:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image4.png|180px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(1)
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Conservation of momentum and mass for incompressible Newtonian fluids in the lagrangian frame of reference are represented by the Navier-Stokes equations and the continuity equation in the final [[Image:Draft_Samper_772971023-image2.png|12px]] position, as follows:
| + | |
| − | | + | |
| − | ''Mass conservation:''
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <math display="inline">\frac{D\rho }{Dt}+\rho \frac{\partial u_i}{\partial x_i}=</math><math>0</math> .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(2)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | ''Momentum conservation:''
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <math display="inline">\rho \frac{Du_i}{Dt}=-\frac{\partial }{\partial x_i}p+</math><math>\frac{\partial }{\partial x_j}{\tau }_{ij}+\rho f_i</math> ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(3)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image7.png|12px]] is the density, [[Image:Draft_Samper_772971023-image8.png|12px]] the pressure, [[Image:Draft_Samper_772971023-image9.png|18px]] the deviatoric stress tensor, [[Image:Draft_Samper_772971023-image10.png|18px]] the source term (usually the gravity) and [[Image:Draft_Samper_772971023-image11.png|24px]] represents the total or material time derivative.
| + | |
| − | | + | |
| − | For Newtonian fluids the stress tensor [[Image:Draft_Samper_772971023-image12.png|18px]] may be expressed as a function of the velocity field through the viscosity [[Image:Draft_Samper_772971023-image13.png|12px]] by
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image14.png|210px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(4)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For near incompressible flows [[Image:Draft_Samper_772971023-image15.png|84px]] the term:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image16.png|78px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(5)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | and it may be neglected from Eq. 4. Then:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image17.png|132px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(6)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In the same way, the term <math display="inline">\frac{\partial }{\partial x_j}{\tau }_{ij}</math> in the momentum equations may be simplified for near incompressible flows as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image19.png|426px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(7)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Using eq. (7), the momentum equations can be finally written as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image20.png|414px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(8)
| + | |
| − | |}
| + | |
| − | | + | |
| − | Note: eq.(3) or the equivalent for incompressible fluid flow eq.(8) are non-linear. In eulerian formulations the non-linearity is explicitly present in the convective terms. In this lagrangian formulation, the non-linearity is due to the fact that eqs. (3) and (8) are written in the final positions of the particles , which are unknown. There are others way to write lagrangian formulations, for instance staying in the initial position [aub04]. In all cases, the equations are non-linear.
| + | |
| − | | + | |
| − | ====Boundary conditions====
| + | |
| − | | + | |
| − | On the boundaries, the standard boundary conditions for the Navier-Stokes equations are:
| + | |
| − | | + | |
| − | {| style="margin: 1em auto 0.1em auto;"
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 67%;"|[[Image:Draft_Samper_772971023-image21.png|108px]]
| + | |
| − | | style="vertical-align: top;"|on [[Image:Draft_Samper_772971023-image22.png|18px]]
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 67%;"|<math display="inline">u_i{\nu }_i={\overline{u}}_n</math>
| + | |
| − | | style="vertical-align: top;"|on [[Image:Draft_Samper_772971023-image24.png|18px]]
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 67%;"|<math display="inline">u_i{\zeta }_i={\overline{u}}_t</math>
| + | |
| − | | style="vertical-align: top;"|on [[Image:Draft_Samper_772971023-image26.png|18px]]
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image27.png|12px]] and [[Image:Draft_Samper_772971023-image28.png|18px]] are the components of the normal and tangent vectors to the boundary.
| + | |
| − | | + | |
| − | ===2.2 Solid dynamics problem: updating the rigid body position===
| + | |
| − | | + | |
| − | In this paper, the structure will be considered as a rigid solid. Then, the equations of motion for a rigid body are:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image29.png|84px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(9)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where ''F<sub>i</sub>'' are the resultant of the external forces (surface forces, gravity force, etc.), whose line of action passes through the mass center of the body, ''U<sub>i</sub>'' is the velocity of the mass center and ''m ''the total mass of the solid.
| + | |
| − | | + | |
| − | The actual motion of the rigid body consists in the superposition of the translation produced by the resultant force ''F<sub>i</sub>'' and the rotation produced by the couple ''T<sub>i</sub>'' satisfying:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image30.png|66px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(10)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where ''M<sub>i</sub>'' is the angular momentum about the mass center. It must be noted that in (10) the time derivative is expressed as the rate of change with respect to any non-rotating system of axis. It may be also expressed as the derivative with respect to the body fixed axes by:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image31.png|216px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(11)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where <math display="inline">\Omega </math> denotes the angular velocity of the body, ''e'' are orthogonal unit basis vectors, <math display="inline">\epsilon </math> the permutation symbol and [[Image:Draft_Samper_772971023-image34.png|48px]] is the derivative with respect to the body fixed axes.
| + | |
| − | | + | |
| − | Let now the body fixed axes be the principal axes of inertia of the body, with its origin at the center of mass, then:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image35.png|66px]] , (without summation in the index ''i'')
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(12)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where ''I<sub>i</sub> ''are the principal moments of inertia and then:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image36.png|108px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(13)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Finally, the equations of motion of the body might be summarized as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image29.png|84px]] ,
| + | |
| − | |}
| + | |
| − | | colspan='2' style="width: 5px;text-align: right;white-space: nowrap;"|(14)
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image37.png|198px]] ,
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(15)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Calling ''a<sub>i</sub>'' and [[Image:Draft_Samper_772971023-image38.png|18px]] the linear and the angular acceleration of the mass center of the body:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image39.png|66px]] ,
| + | |
| − | |}
| + | |
| − | | colspan='2' style="width: 5px;text-align: right;white-space: nowrap;"|(16)
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image40.png|180px]] ,
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(17)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | This is a non-linear system of partial differential equations that has to be linearized for its numerical approximation.
| + | |
| − | | + | |
| − | The final rigid body velocity of an arbitrary point is a combination of both, the linear velocity of the center of mass U<sub>i</sub> and the angular velocity [[Image:Draft_Samper_772971023-image41.png|18px]] according to:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image42.png|132px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(18)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where ''r<sub>i</sub>'' is the distance from the origin of the body axes to an arbitrary point attached to the body. The velocity ''u<sub>i</sub>'' will be used later as a boundary condition for the fluid dynamics problem.
| + | |
| − | | + | |
| − | A very large number of problems involve plane motion. In this case, equation (15) reduces to:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image43.png|114px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(19)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where <math display="inline">\Omega </math> , ''I'' , [[Image:Draft_Samper_772971023-image44.png|12px]] and ''T'' are the planar angular velocity, the moment of inertia, the planar angular acceleration and the external couple respectively.
| + | |
| − | | + | |
| − | ==The discrete fluid dynamics problem ==
| + | |
| − | | + | |
| − | The Navier-Stokes equations present three main difficulties:
| + | |
| − | | + | |
| − | * The equations are time dependent and thus a temporal integration needs to be carried out.
| + | |
| − | | + | |
| − | * A spatial dependency is also present and thus the space will be discretized.
| + | |
| − | | + | |
| − | * Finally, Eq. 3 presents a non-linearity, which must be solved iteratively.
| + | |
| − | | + | |
| − | Each of the above items will be explained and a solution algorithm will be introduced to obtain a final accurate and robust numerical scheme.
| + | |
| − | | + | |
| − | ===3.1 Implicit-explicit time integration===
| + | |
| − | | + | |
| − | Let [[Image:Draft_Samper_772971023-image45.png|12px]] and [[Image:Draft_Samper_772971023-image46.png|24px]] be the initial and final time step. Let [[Image:Draft_Samper_772971023-image47.png|90px]] be the time increment.
| + | |
| − | | + | |
| − | Eq. (8) is integrated implicitly in time as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image48.png|600px]] </big>,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(20)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image49.png|72px]] means [[Image:Draft_Samper_772971023-image50.png|306px]]
| + | |
| − | | + | |
| − | and [[Image:Draft_Samper_772971023-image51.png|90px]] represents the value of the function at time [[Image:Draft_Samper_772971023-image52.png|12px]] but at the final position [[Image:Draft_Samper_772971023-image53.png|12px]] . For simplicity [[Image:Draft_Samper_772971023-image54.png|18px]] will be used instead of [[Image:Draft_Samper_772971023-image55.png|18px]] .
| + | |
| − | | + | |
| − | Only the case of [[Image:Draft_Samper_772971023-image56.png|36px]] (fully implicit scheme) will be considered next. Other values, as for instance [[Image:Draft_Samper_772971023-image57.png|54px]] , may be considered without major changes.
| + | |
| − | | + | |
| − | The time integrated equations become:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image58.png|348px]] .</big>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(21)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The mass conservation is also integrated implicitly by:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image59.png|216px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(22)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The time integration of Eq. (20) presents some difficulties: it is a fully coupled equation involving four degrees of freedom by node. When the fluid is incompressible or nearly incompressible advantages can be taken from the fact that in Eq. (20) the three components of the velocity are only coupled via the pressure. The fractional-step method proposed in [cod01] will be used for the time solution. This basically consists in splitting each time step in two pseudo-time steps. In the first step the implicit part of the pressure is avoided in order to have a decoupled equation in each of the velocity components. The implicit part of the pressure is added at a second step. The fractional-step algorithm for eqs. (21) and (22) is the following:
| + | |
| − | | + | |
| − | ==Split of the momentum equations==
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image60.png|462px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(23)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image61.png|18px]] '' ''are fictitious variables termed fractional velocities defined by the split:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image62.png|288px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(24)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image63.png|216px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(25)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image64.png|90px]] is the value of the pressure at time [[Image:Draft_Samper_772971023-image52.png|12px]] evaluated at the final position and [[Image:Draft_Samper_772971023-image10.png|18px]] is considered constant in time.
| + | |
| − | | + | |
| − | In Eqs. (24) and (25) [[Image:Draft_Samper_772971023-image65.png|12px]] is a parameter giving the amount of pressure splitting, varying between 0 and 1. A larger value of [[Image:Draft_Samper_772971023-image65.png|12px]] means a small pressure split. In this paper [[Image:Draft_Samper_772971023-image65.png|12px]] will be fixed to 0 in order to have the larger pressure split and hence, a better pressure stabilization. Other values as, for instance [[Image:Draft_Samper_772971023-image66.png|36px]] , may be used to derive high order schemes in time[cod01].
| + | |
| − | | + | |
| − | Taking into account Eq. (7), the last term in Eq. (24) may be written as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image67.png|450px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(26)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The following approximations have been introduced [cod01]:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image68.png|360px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(27)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | This allows to write Eq. (24) as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image69.png|498px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(28)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For [[Image:Draft_Samper_772971023-image66.png|36px]] and [[Image:Draft_Samper_772971023-image56.png|36px]] the equations for the fractional velocities becomes:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image70.png|234px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(29)
| + | |
| − | |}
| + | |
| − | | + | |
| − | ==Split of the mass conservation equations==
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image71.png|414px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(30)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image72.png|18px]] is a fictitious variable defined by the split
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | | rowspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|<big> [[Image:Draft_Samper_772971023-image73.png|192px]] </big>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(31)
| + | |
| − | |-
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(32)
| + | |
| − | |}
| + | |
| − | | + | |
| − | ==Coupled equations==
| + | |
| − | | + | |
| − | From eqs. (25), (31) and (32) the coupled mass-momentum equation becomes:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image74.png|168px]] </big>.
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(33)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Taking into account Eq. (31) the above expression can be written as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image75.png|234px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(34)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | It is important to note that in eq. (34) the incompressibility condition has not be introduced yet. The simplest way to introduce the incompressibility condition in a lagrangian formulation is to write:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image76.png|132px]] </big>.
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(35)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Then, the first term of Eq. (34) disappears, giving:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image77.png|150px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(36)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The three steps of the fractional method used here can be summarized by:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image78.png|234px]] [[Image:Draft_Samper_772971023-image79.png|36px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(37)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image77.png|150px]] [[Image:Draft_Samper_772971023-image80.png|48px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(38)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image81.png|174px]] [[Image:Draft_Samper_772971023-image82.png|48px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(39)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | ===3.2 The spatial discretization provided by the MFEM===
| + | |
| − | | + | |
| − | One of the key to solve a fluid mechanics problem using a lagrangian formulation is to generate efficiently the shape functions to approximate the spatial unknown. In the Finite Element context, this means to generate permanently, at each time step, a new mesh. In this work the interpolation function used by the Meshless Finite Element Method (MFEM) [ide03a] will be implemented. This function uses the Voronoï diagram of the cloud of points to construct the interpolant. The extended Delaunay tessellation (EDT) [ide03b] is applied to connect the neighboring particles. The EDT provides polyhedral elements that are sliver-free in 3-D, avoiding instabilities of the Delaunay tessellation due to distorted tetrahedral. EDT provide a way to generate meshes at each time step very efficiently in a computing time which is largely smaller than the computing time needed to solve the linearized system of equation. EDT together with the MFEM are the main key to make the PFEM presented in this paper a useful tool.
| + | |
| − | | + | |
| − | The unknown functions are approximated using an equal order interpolation for all variables in the final configuration:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | | rowspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|<big> [[Image:Draft_Samper_772971023-image83.png|132px]] </big>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(40)
| + | |
| − | |-
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(41)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In matrix form:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image84.png|114px]]
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(42)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image85.png|102px]] ,
| + | |
| − | |}
| + | |
| − | | colspan='2' style="width: 5px;text-align: right;white-space: nowrap;"|(43)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | or in compact form:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image86.png|216px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(44)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where [[Image:Draft_Samper_772971023-image87.png|24px]] are the MFEM shape functions and [[Image:Draft_Samper_772971023-image88.png|30px]] ''' '''the nodal values of the three components of the unknown velocity and the pressure respectively.
| + | |
| − | | + | |
| − | It must be noted that the shape functions [[Image:Draft_Samper_772971023-image89.png|54px]] are functions of the particle coordinates only. Then, the shape functions may change in time following the particle positions.
| + | |
| − | | + | |
| − | During the time step, a mesh update may introduce a change in the shape function definition, which must be taken into account. During the time integration there are two times involved: [[Image:Draft_Samper_772971023-image52.png|12px]] and [[Image:Draft_Samper_772971023-image90.png|24px]] . The following notation will be used to distinguish between [[Image:Draft_Samper_772971023-image91.png|60px]] and [[Image:Draft_Samper_772971023-image92.png|72px]] :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image93.png|96px]] and [[Image:Draft_Samper_772971023-image94.png|120px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(45)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In this work, the following hypothesis will be introduced: There is not mesh update during each time step. This means that if a mesh update is introduced at the beginning of a time step, the same mesh (but deformed) will be kept until the end of the time step.
| + | |
| − | | + | |
| − | Mathematically this means:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image95.png|150px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(46)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Unfortunately, this hypothesis is not always true and this introduces small errors in the computation, which are neglected in this paper.
| + | |
| − | | + | |
| − | Using the Galerkin weighted residual method to solve the splitted equations, the following integrals are obtained:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image96.png|426px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(47)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image97.png|504px]] </big>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(48)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image98.png|474px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(49)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where the boundary conditions have been also spliced and [[Image:Draft_Samper_772971023-image99.png|18px]] is the volume at time [[Image:Draft_Samper_772971023-image100.png|24px]] .
| + | |
| − | | + | |
| − | Integrating by parts some of the terms, the above equations become:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image101.png|402px]]
| + | |
| − | |}
| + | |
| − | | colspan='2' style="width: 5px;text-align: right;white-space: nowrap;"|(50)
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image102.png|438px]]
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(51)
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image103.png|450px]]
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(52)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The essential and natural boundary conditions of Eq. (51) are:
| + | |
| − | | + | |
| − | {| style="width: 36%;margin: 1em auto 0.1em auto;"
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 67%;"|[[Image:Draft_Samper_772971023-image104.png|42px]]
| + | |
| − | | style="vertical-align: top;width: 32%;"|<math display="inline">\mbox{ }on\mbox{ }{\Gamma }_{\sigma }</math>
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;"|<big> [[Image:Draft_Samper_772971023-image106.png|132px]] </big>
| + | |
| − | | style="vertical-align: top;width: 32%;"|<math display="inline">\mbox{ }on\mbox{ }{\Gamma }_u</math>
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where <math display="inline">u_i{}^{n+1}\vert {}_s</math> is the rigid body velocity obtained from Eq. (19).
| + | |
| − | | + | |
| − | ==Discrete equations==
| + | |
| − | | + | |
| − | Using the approximations given be Eqs. (44), (45) and (46) the discrete equations become:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image109.png|480px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(53)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In compact form:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image110.png|330px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(54)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Making use of the approximation described before for [[Image:Draft_Samper_772971023-image111.png|36px]] :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image112.png|450px]] ,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(55)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For [[Image:Draft_Samper_772971023-image113.png|36px]] and [[Image:Draft_Samper_772971023-image114.png|36px]] :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image115.png|222px]] .
| + | |
| − | |}
| + | |
| − | | style="border-left: 1pt solid black;width: 5px;text-align: right;white-space: nowrap;"|(56)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In the same way:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image116.png|516px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(57)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In compact form:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image117.png|234px]] </big>,
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(58)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For [[Image:Draft_Samper_772971023-image113.png|36px]] and [[Image:Draft_Samper_772971023-image114.png|36px]] :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image118.png|174px]] .
| + | |
| − | |}
| + | |
| − | | style="border-left: 1pt solid black;width: 5px;text-align: right;white-space: nowrap;"|(59)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Finally:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image119.png|468px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(60)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | In compact form:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image120.png|246px]] </big>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(61)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For [[Image:Draft_Samper_772971023-image113.png|36px]] and [[Image:Draft_Samper_772971023-image114.png|36px]] :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | <big> [[Image:Draft_Samper_772971023-image121.png|186px]] </big>.
| + | |
| − | |}
| + | |
| − | | style="border-left: 1pt solid black;width: 5px;text-align: right;white-space: nowrap;"|(62)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Where the matrices are:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | | colspan='4' style="vertical-align: top;text-align: center;white-space: nowrap;"|<big> [[Image:Draft_Samper_772971023-image122.png|162px]] </big>,
| + | |
| − | | colspan='4' style="width: 5px;text-align: right;white-space: nowrap;"|(63)
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;text-align: center;white-space: nowrap;"|<big> [[Image:Draft_Samper_772971023-image123.png|114px]] </big>,
| + | |
| − | | colspan='6' style="width: 5px;text-align: right;white-space: nowrap;"|(64)
| + | |
| − | |-
| + | |
| − | | colspan='7' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image124.png|312px]] ,
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(65)
| + | |
| − | |-
| + | |
| − | | colspan='5' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image125.png|270px]] ,
| + | |
| − | | colspan='3' style="width: 5px;text-align: right;white-space: nowrap;"|(66)
| + | |
| − | |-
| + | |
| − | | colspan='3' style="vertical-align: top;text-align: center;white-space: nowrap;"|[[Image:Draft_Samper_772971023-image126.png|120px]] ,
| + | |
| − | | colspan='5' style="width: 5px;text-align: right;white-space: nowrap;"|(67)
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image127.png|108px]] ,
| + | |
| − | |}
| + | |
| − | | colspan='7' style="width: 5px;text-align: right;white-space: nowrap;"|(68)
| + | |
| − | |-
| + | |
| − | | colspan='6' style="vertical-align: top;text-align: center;white-space: nowrap;"|<big> [[Image:Draft_Samper_772971023-image128.png|324px]] </big>
| + | |
| − | | colspan='2' style="width: 5px;text-align: right;white-space: nowrap;"|(69)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | ====Stabilization of the incompressibility condition====
| + | |
| − | | + | |
| − | In the eulerian form of the momentum equations, the discrete form must be stabilized in order to avoid numerical wiggles in the velocity and pressure results. This is not the case in the lagrangian formulation where no stabilization terms must be added to eqs. (62). Nevertheless, the incompressibility condition must be stabilized in equal-order approximations to avoid pressure oscillations in some particular cases.
| + | |
| − | | + | |
| − | For instance for small pressure split ( [[Image:Draft_Samper_772971023-image129.png|36px]] ) or for small time step increments (Courant number much less than one) it is well known that the fractional step does not stabilize the pressure waves [cod01]. In those particular cases, a stabilization term must be introduced in Eqs. (62) in order to eliminate pressure oscillations.
| + | |
| − | | + | |
| − | A simple and effective procedure to derive a stabilized formulation for incompressible flows is based in the so-called Finite Calculus (FIC) formulations [ona98, ona00, ona02].
| + | |
| − | | + | |
| − | In all the examples presented in this paper the FIC formulation was used to stabilize the pressure oscillations.
| + | |
| − | | + | |
| − | ===3.3 Non-linearity of the lagrangian formulation===
| + | |
| − | | + | |
| − | Many algorithms are available to linearize the equations of motion. The Newton-Raphson scheme is probably the most popular because of its robustness and fast convergence. It has been applied with success in this type of formulation in [rad98]. Nevertheless we consider that it might not be the most appropriate option for the type of equations we are intending to solve as it requires large memory storage. Instead, the successive iteration algorithm has been chosen for the present analysis. In this case, only the variables that induces the non-linearity need to be stored in successive iterations. Let us now describe the process that may take place until convergence:
| + | |
| − | | + | |
| − | :1) Approximate [[Image:Draft_Samper_772971023-image130.png|24px]] (For the first iteration [[Image:Draft_Samper_772971023-image131.png|54px]] . For the subsequent iterations the value of [[Image:Draft_Samper_772971023-image130.png|24px]] corresponding to the last iteration is taken).
| + | |
| − | | + | |
| − | :2) Move the particles to the [[Image:Draft_Samper_772971023-image132.png|24px]] position and perform an EDT polyhedrization.
| + | |
| − | | + | |
| − | :3) Evaluate the fractional velocity [[Image:Draft_Samper_772971023-image133.png|18px]] from (56). It must be noted that the matrices [[Image:Draft_Samper_772971023-image134.png|18px]] and [[Image:Draft_Samper_772971023-image135.png|18px]] are separated in 3 blocks. Then, these equations may be solved separately for [[Image:Draft_Samper_772971023-image136.png|48px]] and [[Image:Draft_Samper_772971023-image137.png|18px]] . For [[Image:Draft_Samper_772971023-image138.png|36px]] (implicit scheme) involves the solution of 3 symmetric linear systems of equations. For [[Image:Draft_Samper_772971023-image139.png|36px]] (explicit scheme) the M matrix may be lumped and inverted directly.
| + | |
| − | | + | |
| − | :4) Evaluate the pressure [[Image:Draft_Samper_772971023-image140.png|30px]] by solving the laplacian equation (58).
| + | |
| − | | + | |
| − | :<big>5) </big>Evaluate [[Image:Draft_Samper_772971023-image141.png|24px]] using (61).
| + | |
| − | | + | |
| − | == Time integration of the solid dynamics problem==
| + | |
| − | | + | |
| − | Eqs. (14) and (15) that govern the movement of rigid bodies are integrated in time by the explicit Newmark algorithm. It consists in evaluating the velocity by linearizing the acceleration between two time steps:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image142.png|222px]] .
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image143.png|228px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(70)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The point position for the explicit version of the Newmark algorithm is evaluated by:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image144.png|198px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(71)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | To integrate the angular acceleration in a 3D system by Newmark algorithm two steps are needed, namely:
| + | |
| − | | + | |
| − | ==predictor step:==
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image145.png|162px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(72)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | then, the accelerations are predicted by using Eq. (17):
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image146.png|240px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(73)
| + | |
| − | |}
| + | |
| − | | + | |
| − | ==corrector step==
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image147.png|144px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(74)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The linear velocities are integrated directly using Eqs.(16) and (70):
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image148.png|102px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image142.png|222px]] .
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(75)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | For the present analysis [[Image:Draft_Samper_772971023-image149.png|54px]] will be considered.
| + | |
| − | | + | |
| − | Both velocities, [[Image:Draft_Samper_772971023-image150.png|30px]] and [[Image:Draft_Samper_772971023-image151.png|30px]] are used in (18) to evaluate the new velocity of all the points of the body.
| + | |
| − | | + | |
| − | In the explicit version of the Newmark algorithm, the new position of the rigid body is evaluated by :
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image144.png|198px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image152.png|168px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;"|(76)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | where ''x<sub>i </sub>''is the new position of the center of mass and [[Image:Draft_Samper_772971023-image153.png|12px]] <sub>i</sub> are the angular rotations of the body.
| + | |
| − | | + | |
| − | It must be noted that for planar motion, the predictor step is unnecessary and [[Image:Draft_Samper_772971023-image154.png|30px]] may be evaluated directly using (19).
| + | |
| − | | + | |
| − | ===4.1 The Coupled Problem===
| + | |
| − | | + | |
| − | On the coupling boundary the fluid velocity and the solid velocity should converge to the same value. This could be expressed as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image155.png|84px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (77)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | Thus, two subsystems need to be solved, namely the fluid system:
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image156.png|210px]]
| + | |
| − | | + | |
| − | and the solid system:
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image157.png|204px]] .
| + | |
| − | | + | |
| − | In the equations above only the variables to be solved at time step [[Image:Draft_Samper_772971023-image158.png|30px]] are shown.
| + | |
| − | | + | |
| − | An iterative procedure has to be implemented to couple both systems. A fixed point algorithm may be implemented and thus the system could be written as:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: center;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: center; margin:auto;"
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image159.png|312px]]
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (78)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | The first equation in Eqs. (78) denotes the fluid subsystem and the second equation the solid subsystem. The subscript ''k'' is the iteration counter.
| + | |
| − | | + | |
| − | In the present analysis a Gauss-Seidel process has been chosen [cod96]. In this way, the iterative procedure means to solve first one of the subsystems, for instance the fluid system. Next the solid system is solved using the information from the fluid computation. Eqs. (78) should be modified and the final expression used for the computation is as follow:
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image160.png|312px]]
| + | |
| − | | + | |
| − | Convergence occurs when the difference between velocities of successive iteration steps is less than the acceptable error.
| + | |
| − | | + | |
| − | == Free-surface and boundary recognition==
| + | |
| − | | + | |
| − | The solution of partial differential equations (PDE) requires to prescribe boundary conditions as a necessary step to a well-posed problem. When the PDEs are approximated in space and the domain is partitioned into discrete elements (finite elements, particles, balls, nodes, etc.) the boundary elements should be provided at the initial time step, such that, at run time the algorithm knows where to impose or fix the variables of the analysis (pressure, velocity and their derivatives for instance). This would be the case of a static domain, where the geometry does not change in time and the boundaries remain constant.
| + | |
| − | | + | |
| − | In this work the interest is focused on problems where the solution domain is highly distorted, and boundary elements can change between time steps. In this case an efficient boundary recognition algorithm is mandatory in order to impose boundary conditions over the right elements, thus avoiding possible error accumulation over the time.
| + | |
| − | | + | |
| − | When applying the MFEM [Ide03a] to the discrete space problem, the EDT [Ide03b] is computed to connect the particles that discretize the domain, thus, all the empty Voronoï spheres are found and stored. These spheres will be used to compute the boundary using the alpha-shape technique [ede94].
| + | |
| − | | + | |
| − | The particles will follow a given ''h(x)'' distribution according to the maximum error allowed for the discrete space problem, where ''h(x)'' is the expected distance among neighboring particles. Then, having all the empty Voronoï spheres and ''h(x)'' the boundary particles are regarded as: ''all the particles which are on an empty sphere with a radius [[Image:Draft_Samper_772971023-image161.png|12px]] bigger than '' [[Image:Draft_Samper_772971023-image162.png|24px]] .
| + | |
| − | | + | |
| − | In this criterion, [[Image:Draft_Samper_772971023-image163.png|12px]] is a parameter close to one, typically [[Image:Draft_Samper_772971023-image164.png|48px]] and [[Image:Draft_Samper_772971023-image165.png|12px]] is the mean value taken from the defining particles of the sphere under inspection.
| + | |
| − | | + | |
| − | Once a decision has been made concerning which of the nodes are on the boundaries, the boundary surface must be defined. It is well known that, in 3-D problems, the surface fitting a number of nodes is not unique. For instance, four boundary nodes on the same sphere may define two different boundary surfaces, one concave and the other convex.
| + | |
| − | | + | |
| − | In this work, the boundary surface is defined by all the polyhedron faces having all their nodes on the boundary and belonging to just one polyhedron. See [ide03b].
| + | |
| − | | + | |
| − | The correct boundary surface may be important to define the correct normal external to the surface. Furthermore; in weak forms (Galerkin) it is also important a correct evaluation of the volume domain. It must however be noted that in the criterion proposed above, the error in the boundary surface definition is of order [[Image:Draft_Samper_772971023-image166.png|12px]] . This is the standard error of the boundary surface definition in a meshless method for a given node distribution.
| + | |
| − | | + | |
| − | Another important feature of the alpha-shape technique related to the contact problem is shown in Fig.5.1. The image shows two different time steps of a solid cube falling into water after the alpha-shape algorithm for the boundary recognition has been applied. The particles, as well as the connections provided by the EDT are depicted. At the first time step all the radii of the empty circles constructed with the nodes of the cube and the nodes regarded as belonging to the free-surface are larger than [[Image:Draft_Samper_772971023-image167.png|48px]] and thus the elements that they define are eliminated from the tessellation. The second picture shows a more evolved state, with the cube reaching the water surface. Thus, at this state the circle radii are less than [[Image:Draft_Samper_772971023-image167.png|48px]] and so the connections between the cube and the fluid take part of the computation. In this way free surface and contacts are solved at once.
| + | |
| − | | + | |
| − | {| style="width: 100%;margin: 1em auto 0.1em auto;"
| + | |
| − | |-
| + | |
| − | | style="vertical-align: top;width: 50%;"|[[Image:Draft_Samper_772971023-image168.jpeg|198px]]
| + | |
| − | | style="vertical-align: top;width: 49%;"|[[Image:Draft_Samper_772971023-image169.jpeg|198px]]
| + | |
| − | |-
| + | |
| − | | colspan='2' style="vertical-align: top;"|'''Fig. 5.1''' '': The alpha-shape technique used for contact recognition.''
| + | |
| − | |}
| + | |
| − | | + | |
| − | =6. Joining and breaking particles=
| + | |
| − | | + | |
| − | The idea of ''h'' variable mesh is rather different in particle methods than in classical eulerian formulations. In particle methods, each particle is followed in time and the same particle can cross domains in which the solution need small ''h'' in order to represent high gradients or can cross a region with large ''h'' where the solution is smooth. The concept of variable ''h'' is introduce in particle methods by joining two particles when they are too close to each other or breaking a particle in two when all the neighboring particles are too far and the solution needs a higher gradient.
| + | |
| − | | + | |
| − | In the example presented in this paper the following criterion has been used:
| + | |
| − | | + | |
| − | 1) During the EDT algorithm to build the polyhedral mesh a particle is not added if there is a previous point at a distance ''d'' < λ ''h''(''x''), being λ = 0.5 a constant parameter.
| + | |
| − | | + | |
| − | 2) on the contrary, when there is an empty sphere whose radius r > ''h''(''x'') a point is added to its center and assigned with values interpolated from the sphere defining particles. The parameter is currently taken as = 1.1*sqrt(dim)/2 in order to accept, with a 10% of tolerance, near-square (dim = 2, ≈ .78) or near-cube (dim = 3, ≈ .95) local arrays as connected inner points. This parameter must be less than the alpha-shape parameter [[Image:Draft_Samper_772971023-image167.png|48px]] in order to avoid interference with the boundary recognition algorithm.
| + | |
| − | | + | |
| − | ==7. Validation Examples.==
| + | |
| − | | + | |
| − | PFEM was developed as a general-purpose method for solving different kind of problems on which large free surface or interface boundaries changes are involved. The method is well suited to solve a large variety of mechanical problems including mixing fluid and solid materials, wave motion problems, mould filling, coupled thermal-mechanical problems and fluid-solid interaction as well.
| + | |
| − | | + | |
| − | In this section some problems are included in order to show the validation of the present approach towards experimental and numerical tests. In the following section, some more specific examples on fluid structure interaction will be shown.
| + | |
| − | | + | |
| − | ===7.1 Sloshing===
| + | |
| − | | + | |
| − | The simple problem of the free oscillation of an incompressible liquid in a container is considered first. Numerical solutions for this problem can be found in several references [rad98]. This problem is interesting because there is an analytical solution for small amplitudes. For larger amplitudes the wave breaks and also some particles can be separated from the fluid domain due to their large velocity. PFEM can solve very large amplitudes, even in a 3D domain [ide04]. However, in this section only the small amplitude, twodimensional example is shown to validate the method.
| + | |
| − | | + | |
| − | Figure 7.1 shows the variation in time of the amplitude compared with the analytical results for the near inviscid case. Little numerical viscosity is observed on the phase wave and amplitude in spite of the relative poor point distribution.
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image170.png|486px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 7.1 Sloshing: Comparison of the numerical and analytical solution.'''</div>
| + | |
| − | | + | |
| − | ===7.2 Wave on a channel: comparison with experimental results.===
| + | |
| − | | + | |
| − | This example was performed in order to compare and validate the method with experimental results. A wave is running from the left to the right arriving to a shallow domain were the wave breaks. The example is represented on Fig 7.2 were the calculated particle positions are shown at different time steps. The wave was produce by a particular movement of the left wall. This example was reproduced experimentally in the CIEM (Maritime Experimental and Research Channel) of the Escuela Técnica de Ingenieros de Canales Caminos y Puertos in the University of Catalunya. The channel is 100 meters length, 3 meters wide and 5 meters high.
| + | |
| − | | + | |
| − | A pressure sensor was placed on the right wall at 0.2 meters from the bottom.
| + | |
| − | | + | |
| − | <br/> [[Image:Draft_Samper_772971023-image171-c.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image172.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image173.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image174-c.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image175-c.png|600px]] [[Image:Draft_Samper_772971023-image176-c.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image177-c.png|600px]]
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 7.2''' '' Wave on a channel: particle distribution for different time steps''.</div>
| + | |
| − | | + | |
| − | Experimental and numerical pressure results are compared at different time step in Fig. 7.3. Both results, experimental and numerical, were smoothed in order to ignore the artificial oscillations from high order waves present in the problem. The comparison of the results in the pressure values shows a reasonably agreement
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-chart1.svg|600px]]
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 7.3''' '' Wave on a channel: experimental and numerical comparison of the pressure''.</div>
| + | |
| − | | + | |
| − | ===7.3 Dam collapse===
| + | |
| − | | + | |
| − | | + | |
| − | {|
| + | |
| − | |-
| + | |
| − | | [[Image:Draft_Samper_772971023-image178.png|264px]]
| + | |
| − | | [[Image:Draft_Samper_772971023-image179.png|center|252px]]
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 7.4.''' '' Dam Collapse. Initial position. Left: experimental [Kos96]. Right: 3D simulation.''</div>
| + | |
| − | | + | |
| − | The dam collapse problem represented in Figure 7.4 was solved by Koshizu and Oka [Kos96] both experimentally and numerically in a 2D domain. It became a classical example to test the validation of the Lagrangian formulation in fluid flows.
| + | |
| − | | + | |
| − | The water is initially located on the left supported by a removable board. The collapse starts at time ''t ''= 0, when the removable board is slid-up. Viscosity and surface tension are neglected. The water is running on the bottom wall until, near 0.3 sec, it impinges on the right vertical wall. Breaking waves appear at 0.6 sec. Around t=1 sec. the main water wave reaches the left wall again.
| + | |
| − | | + | |
| − | In [Ide04] the results obtained using the method proposed in 2D and 3D domains are presented and compared with experimental results. Agreement with the experimental results of [kos96] both in the shape of the free surface and in the time development are excellent.
| + | |
| − | | + | |
| − | In this example the power of the method to represent breaking waves and flow separation for a very complicated and random problem is verified and compared to experimental results.
| + | |
| − | | + | |
| − | This example is further exploited here to compare results obtained with three different node densities at some time steps in order to check the convergence of the method.
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image180.png|600px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 7.5.''' '' Dam Collapse. Comparison between the experiment and numerical results obtained in different time steps, with different refinement levels.''</div>
| + | |
| − | | + | |
| − | Figure 7.5 shows the domain profile at different time steps and with different refinement levels. At the top there are some photos taken on the experimental setup corresponding to .2, .6 and 1 second from the starting point. From top to bottom grids with 2.2, 4.4 and 8.8mm in typical distance between neighboring particles are shown.
| + | |
| − | | + | |
| − | The excellent agreement already obtained with low refinement is even better with higher level of refinement.
| + | |
| − | | + | |
| − | ==8 Further Examples on fluidstructure interaction==
| + | |
| − | | + | |
| − | ===8.1 Ship profile hit by a wave===
| + | |
| − | | + | |
| − | In the example of Figure 8.1 the motion of a fictitious rigid ship hit by an incoming wave is analyzed. This is the first example in which the rigid body is moved by the fluid forces in a coupling problem as was explained in the previous chapters.
| + | |
| − | | + | |
| − | The dynamic motion of the ship is induced by the resultant of the pressure and the viscous forces acting on the ship boundaries. The horizontal displacement of the mass centre of the ship was fixed to zero. In this way, the ship moves only vertically although it can freely rotate. The position of the ship boundary at each time step is evaluated using eq. (76) and the velocity of the body by using eq. (18). This defines a moving boundary condition for the free surface particles in the fluid as introduced in eq. (59).
| + | |
| − | | + | |
| − | Figure 7.7 shows instants of the motion of the ship and the surrounding fluid. It is interesting to see the breaking of a wave on the ship prow at t= 0.91 sec. as well as on the stern at t= 2.05 sec. when the wave goes back. Note that some water drops slip over the ship deck at t= 1.3 sec. and 2.95 sec.
| + | |
| − | | + | |
| − | ==8.2 Oil ship tank under a lateral wave.==
| + | |
| − | | + | |
| − | The present example depicts the flexibility of the algorithm introduced in this paper to solve some complicated configurations as the one shown in Fig. 8.2. The traversal cut of an oil ship tank is hit by a wave. The structure of the ship does not only interact with the external water but it also moves due to the fluid forces induced by the fluid in the tank.
| + | |
| − | | + | |
| − | Fig. 8.2 shows the temporal development of the problem. The blue lines over each particle represent the magnitude of the velocity field.
| + | |
| − | | + | |
| − | Initially (t = 0.0) the ship is released from a fixed position and the equilibrium configuration found is consistent with Arquimides principle. During the following time steps the external wave hits the ship and both the internal and the external fluids interact with the ship boundaries. At times t = 5.10 sec and 6.00 sec. breaking waves and some water drops slipping along the ship deck can be observed. Figure 8.3 shows the pressure pattern at two time steps.
| + | |
| − | | + | |
| − | ===8.3. Tanker Sinking===
| + | |
| − | | + | |
| − | This example represents the sinking of a tanker, which is being flooded by the ship prow. The ship has many connected compartments that are serially flooded.
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | <br/> [[Image:Draft_Samper_772971023-image181.png|600px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image182.png|600px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image183.png|600px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image184.png|600px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image185.png|516px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image186.png|516px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 8.1.''' '' Ship profile hit by a wave: particle positions for different time steps''</div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image187-c.png|486px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image188-c.png|486px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 8.2:''' '' Oil ship tank under a lateral wave: Particle distribution and velocity field''.</div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image189-c.png|492px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | <big> [[Image:Draft_Samper_772971023-image190-c.png|492px]] </big></div>
| + | |
| − | | + | |
| − | '''Fig. 8.3:''' '' Oil ship tank under a lateral wave: Pressure profile at two different time steps.''
| + | |
| − | | + | |
| − | <br/> [[Image:Draft_Samper_772971023-image191.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image192.png|600px]]
| + | |
| − | | + | |
| − | [[Image:Draft_Samper_772971023-image193.png|600px]]
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 8.4:''' '' Sinking tanker. Particle distribution at three different timesteps.''</div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | <span style="text-align: center; font-size: 75%;"> [[Image:Draft_Samper_772971023-image194.png|456px]] </span></div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | [[Image:Draft_Samper_772971023-image195.png|438px]] </div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | <span style="text-align: center; font-size: 75%;"> [[Image:Draft_Samper_772971023-image196.png|498px]] </span></div>
| + | |
| − | | + | |
| − | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
| + | |
| − | '''Fig. 8.5:''' '' Sinking tanker. Velocity field.''</div>
| + | |
| − | | + | |
| − | In this example proper fluid-structure interaction is displayed as the buoyancy, pressure and drag forces from the fluid are acting over the ship and, on the other hand, the ship displacement moves internal and external free surfaces on the fluid.
| + | |
| − | | + | |
| − | In Figure 8.4 there are three timesteps shown with the particle positions and the tanker in three different sinking stages.
| + | |
| − | | + | |
| − | Figure 8.5 displays velocity profile on a zoom of the first and last time steps form previous figure. In this figure, the vorticity is also easily shown.
| + | |
| − | | + | |
| − | This is an interesting example using a variable distance between particles to enhance the solution near the ship and free surfaces. This variable distribution was obtained following the method outlined in section 6 above.
| + | |
| − | | + | |
| − | The large cyan dots are representing the free surface recognized by the method as explained in section 5.
| + | |
| − | | + | |
| − | ==9. Conclusions==
| + | |
| | | | |
| − | The particle finite element method (PFEM) seems ideal to treat problems involving fluids with free-surface and submerged or floating structures within a unified lagrangian finite element framework. Problems such as the analysis of fluid-structure interactions, large motion of fluid or solid particles, surface waves, water splashing, separation of water drops, etc. can be easily solved with the PFEM.
| + | '''Keywords: '''Fluid-Structure interaction, Particle methods, Lagrange formulations, Incompressible Fluid Flows, Meshless Methods, Finite Element Method. |
| | | | |
| − | The success of the method lays in the accurate and efficient solution at each time step of the equations of an incompressible fluid and the interacting solid. Essential ingredients of the numerical solution are the efficient regeneration of the polyhedral mesh using an extended Delaunay tessellation, the polyhedral finite element interpolation via the MFEM and the identification of the boundary nodes using an Alpha-Shape type technique.
| |
| | | | |
| − | The examples presented have shown the potential of the PFEM for solving a wide class of practical FSI problems. Other examples of application of the PFEM can be found in [ona04].
| + | <pdf>Media:Draft_Samper_772971023_6894_con230bis.pdf</pdf> |
| | | | |
| − | ==Acknowledgements==
| |
| | | | |
| − | Thanks are given to Mr. Miguel Angel Calena for computing the numerical results using the PFEM of the wave on a channel presented in 7.2 Thanks also to the Maritime Experimental and Research Channel (CIEM) of the Escuela Técnica de Ingenieros de Canales Caminos y Puertos, University of Catalunya for the experimental results of the same problem.
| |
| | | | |
| | ==References== | | ==References== |
In the present work a new approach to solve fluid-structure interaction problems is described. Both, the equations of motion for fluids and for solids have been approximated using a material (lagrangian) formulation. To approximate the partial differential equations representing the fluid motion, the shape functions introduced by the Meshless Finite Element Method (MFEM) have been used. Thus, the continuum is discretized into particles that move under body forces (gravity) and surface forces (due to the interaction with neighboring particles). All the physical properties such as density, viscosity, conductivity, etc., as well as the variables that define the temporal state such as velocity and position and also other variables like temperature are assigned to the particles and are transported with the particle motion. The so called Particle Finite Element Method (PFEM) provides a very advantageous and efficient way for solving contact and free-surface problems, highly simplifying the treatment of fluid-structure interactions.
[aub04] Aubry R., Idelsohn S.R. and Oñate. “The Particle Finite Element Method in
[bel94] T. Belytschko, L. Gu and Y. Y. Lu, “Fracture and crack growth by element free Galergin methods”, Modelling Simul. Mater Sci. Eng., vol. 2, pp. 519-534, 1994.
[cod96] Codina R, Cervera M. “Block-iterative algorithms for nonlinear coupled problems”. In: Papadrakakis M, Bugeda G, editors, Advanced Computational Methods in Structural mechanics. Barcelona: CIMNE, 1996.
[dua95] C. Duarte and J. D. Oden, “Hp clouds-a meshless method to solve boundary value problems”, Tech. Rep. 95-05, Texas institute for Computational and Applied Mechanics, Univerfsity of Texas, 1995.
[ede94] H. Edelsbrunner and E.P. Mucke, “Three-dimensional alpha-shape”. ACM Transactions on Graphics, 3, (1994), 43–72.
[ide03a] S.R. Idelsohn, E. Oñate, N. Calvo and F. Del Pin, “The meshless finite element method”, Int. J. for Numerical Methods in Engineering, Vol 58, Issue 4, (2003).
[ide04] S.R. Idelsohn, E. Oñate and F. Del Pin, “The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves”. Int. J. for Numerical Methods in Engineering, Vol 61, Issue 7, (2004), Pages: 964-989.
[ide03b] S.R. Idelsohn, N. Calvo and E. Oñate, “Polyhedrization of an arbitrary 3D point set”. Computer Method in Applied Mechanics and Engineering, 192, pp 2649-2667 (2003).
[kos96] S. Koshizuka and Y. Oka “Moving particle semi-implicit method for fragmentation of incompressible fluid”, Nuclear Engineering Science, 123, (1996), 421–434.
[loh03] R. Löhner, J.D. Baum, E.L. Mestreau, D. Sharov, Ch. Charman and D. Pelessone “Adaptive Embedded Unstructured Grid Methods”; AIAA-03-1116 (2003).
[mok01] Mok DP, Wall WA. “Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures”. In: Wall WA, Bletzinger K-U, Schweizerhof K, editors. Trends in Computational Structural Mechanics. Barcelona: CIMNE, 2001.
[mon81] R.A. Gingold and J.J. Monaghan, “Kernel estimates as a basis for general particle methods in hydrodynamics”, Journal of Computational Physics, vol. 46, pp. 429-453, 1981.
[mon97] R.A. Gingold and J.J. Monaghan, “Smoothed particle hydrodynamics, theory and application to non-spherical stars”, Mon. Nat. Roy. Astr. Soc., 181, (1997), 375–389.
[nay92] B. Nayroles, G. Touzot and P. Villon. “Generalizing the finite element method: diffuse approximation and diffuse elements”, Computational Mechanics, vol. 10, pp. 307-318, 1992.
[ona98] E. Oñate. “Derivation of stabilized equations for advective-diffusive transport and fluid flow problems”. Comput. Methods Appl. Mech. Engrg.,1998, 151:1-2, 233--267
[ona00] E. Oñate. “A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation”. Comp. Meth. Appl. Mech. Engng. ,2000, 182:1-2, 355-370.
[ona02] E. Oñate, “Possibilities of finite calculus in computational mechanics”. (2002) Submitted to Int. J. Num. Meth. Engng.
[ona04] E. Oñate, S.R. Idelsohn, F. Del Pin and R. Aubry. “ The Particle Finite Element Method: an overview”. In press in International Journal of Computational Methods. (2004)
[rad98] R. Radovitzky, M. Ortiz “Lagrangian finite element analysis of a Newtonian flows”. Int. J. Numer. Engng. 43, 607-619 (1998).
[ram87] B. Ramaswamy and M. Kawahara, “Lagrangian finite element analysis of Newtonian fluid flows”, Int. J. Num. Meth. Fluids, vol. 7, pp. 953-984, 1987.
[rug00] Rugonyi S, Bathe KJ. “On the analysis of fully-coupled fluid flows with structural interactions - a Coupling and Condensation Procedure”. Int. J of Comp Civil and Struct Engineering 2000; 1:29-41.
[rug01] Rugonyi S, Bathe KJ. “On finite element analysis of fluid flows fully coupled with structural interactions”. Computer Modeling and Simulation in Engineering (CMES) 2001; 2:195-212.
[Sou00] M. Souli, A. Ouahsine and L. Lewin, “Arbitrary Lagrangian Eulerian formulation for Fluid-Structure Interaction Problems”. Comp. Meth. In Applied Mech. and Eng. Vol. 190, pp 659-675 (2000).
[suk98] N. Sukumar, B. Moran and T. Belytschko, “The natural element method in solid mechanics”, International Journal for Numerical Methods in Engineering, vol. 43, pp. 839-887, 1998.
[zha01] Zhang H, Bathe KJ. “Direct and iterative computing of fluid flows fully coupled with structures”, in Computational Fluid and Solid Mechanics (K.J. Bathe, ed.), Elsevier, 2001.
