Benoit Fabrèges

Presentation of the computational code

This code, that I developed during my PhD in collaboration with Loïc Gouarin and Bertrand Maury, is written in C++ and uses the library PETSc. It works with the 4Q1/Q1 ( finite element to simulate the motion of immersed rigid bodies with the smooth extension method developed in my PhD. The library PETSc is used for the parallel structure of matrices and vectors, aswell as the parallel solvers for the linear systems. The code is also interfaced with the code SCoPI in order to handle contacts between the rigid bodies.

Algorithm

The idea of the method is to seek a smooth extension of the exact solution in the whole fictitious domain. For this purpose, the extension of the right-hand side within the rigid bodies is used. The problem amounts to look for an extension of the right-hand side such that the restriction of the solution of a classical Stokes problem, defined in the whole fictitious domain, to the fluid domain is the exact solution of the original problem (the fluid velocity satisfies the rigid motion constraint on the boundary of the rigid bodies). To find such an extension of the right-hand side, a cost function is minimized with a conjugate gradient like algorithm. This minimization process involves the resolution of two classical Stokes problem at each iteration of the conjugate gradient. The algorithm of the code interfaced with the contact code SCoPI is the following:

  • SCoPI sends the position of the rigid bodies to the code.
  • Resolution of a classical Stokes problem to initialize the conjugate gradient.
  • Conjugate gradient loop:
    • Resolution of two classical Stokes problem to get the gradient of the cost function.
    • Update of the extension of the right-hand side
  • End of conjugate gradient
  • Resolution of a classical Stokes problem with the extension of the right-hand side found by the minimization process.
  • Computation of the velocity of each rigid body from the fluid velocity
  • The code sends these velocities to SCoPI that projects them on the space of admissible velocities.

A lot of Stokes problem have to be solved with this method but the only part that is changing from one Stokes problem to another is the right-hand side. The matrices are the one of a classical incompressible Stokes problem and are always the same for all them. Moreover, because we are using cartesian meshes, fast iterative solvers like geometric multigrids are used to solve these Stokes problems.

Some simulations

7500 rigid bodies sedimenting in 3D. The rigid bodies are evenly distributed at initial time and have a radius of 0.0125. The cartesian meshes has the size 256x256x256 for the velocity. The simulation was ran on 24 cores and took 22 minutes by time iteration.
2000 rigid bodies sedimenting in 3D. The top of the cube is rotating. The rigid bodies have different radius. The cartesian mesh has the size 65x65x65 for the velocity.
Sedimentation of 1008 rigid bodies with the same radius of 0.01. The cartesian mesh has a size of 1024x1024 for the velocity.

Additional information

The library FELiScE

The library FELiScE is a finite element library developed at Inria in the teams REO (Paris-Rocquencourt) and MACS (Saclay). This is a parallel code (MPI) written in C++ that uses the library PETSc. In the REO team, we developed it to solve, in particular, fluid dynamics problems, electrophysiology problems or fluid-structure interaction problems. Different methods are implemented in FELiScE for fluid-structure interactions. The one we are interested in here is the Nitsche-XFEM method that we use to simulate the interaction between a fluid and an immersed interface.

The Nitsche-XFEM method

The Nitsche-XFEM method is a fictitious domain method that is able to well-capture strong discontinuities (pressure) or weak ones (fluid velocity) while not suffering from loss of accuracy. Thanks to the fictitious domain method, the fluid and interface meshes are not fitted and are fixed. Large displacement can be simulated without having to deal with deformations or remeshing of the fluid mesh. The fictitious domain is given by the XFEM part of the method that duplicates the fluid mesh elements intersected by the interface. The coupling relation between the fluid and solid equations are enforced accurately with a Nitsche's formulation.

Thanks to this method, we can use P1/P1 finite element for the fluid problem and still capture discontinuities across the interface in both the pressure and the derivatives of the fluid velocity. The implementation of the Nitsche-XFEM method in FELiScE is interfaced with a mesh intersection algorithm developed by Frédéric Alauzet that gives all the information about the intersection between the fluid fixed mesh and the solid mesh. In particular, it is used in the two main key implementation points of the method:

  • Integration over cut-elements. The Nitsche's formulation requires integration in the physical domain only, which is done by integrating a sub-triangulation of each intersected element. One should note that this mesh intersection is not a remeshing of the fluid mesh but only a sub-triangulation used to compute these integrals in the physical domain. The mesh resulting from the mesh intersection algorithm contains some very small and stretched elements that are not suited to run finite element computation but that are perfectly fine to only compute integrals.
Mesh intersection
  • The XFEM method. The degrees of freedom of the fluid mesh elements intersected by the interface are duplicated in order to be able to capture discontinuities. The implementation in FELiScE is general enough to manage the case of an interface not cutting entirely the fluid mesh. In this case, a special treatment is done in the elements containing the tips of the interface. The solution defined in the whole domain (physical + fictitious) is saved so that one can see the smooth extension of the exact solution within the fictitious domain computed with the Nitsche's formulation. In particular, one can note the pressure discontinuities aswell as the discontinuities in the fluid velocity derivatives. It is clear with the velocity elevations that the discontinuities are not fitted to the fluid mesh.
Pressure elevation

Pressure elevation

Pressure elevation

Pressure elevation

The pressure is double valued in the triangle intersected by the interfaces.

Horizontal velocity elevation

Horizontal velocity elevation

Horizontal velocity elevation

Horizontal velocity elevation

Note the jump in the derivatives across the interfaces. The velocity in the physical domain is smoothly extended to the fictitious domain.

Vertical velocity elevation

Vertical velocity elevation

Vertical velocity elevation

Vertical velocity elevation

Note the jump in the derivatives across the interfaces. The velocity in the physical domain is smoothly extended to the fictitious domain.

Our implementation is also flexible. Solid models can be plugged to the Nitsche-XFEM methods easily. Indeed, the class implementing the solid problem in the Nitsche-XFEM method is a kind of decorator whose only aim is to add the specific terms coming from the Nitsche-XFEM formulation. The solid model does not need to be known.

Some videos

A double pressure-wave (fluid from each side of the interface). A sine pressure is imposed on only one side of the interface and then vanishes to let the pressure wave evolve. An elevation of the pressure is represented. One can note the discontinuity in the pressure. A simple elastic string model is used for the interface.
Here, the interface is not cutting the fluid domain in two. There are two valves with no contact and the velocity of the fluid is imposed on the left boundary. Up is the fluid velocity field and bottom is an elevation of the pressure. A linear curved beam model is used for the interfaces.
This is a preliminary result, the parameter are not reflecting any physical experiment. A bubble is entering a canal and deforms. The solid model is again a linear curved beam model and the fluid inside the bubble is the same as the one outside.

Additional informations

One can have more details about the Nitsche-XFEM method in the following submitted article:

Under construction