# Mathematical and computational framework for moving and colliding rigid bodies in a Newtonian fluid

Céline Van Landeghem, Luca Berti, Vincent Chabannes, Christophe

Prud'Homme, Agathe Chouippe, Yannick Hoarau, Laëtitia Giraldi

► To cite this version:

Céline Van Landeghem, Luca Berti, Vincent Chabannes, Christophe Prud'Homme, Agathe Chouippe, et al.. Mathematical and computational framework for moving and colliding rigid bodies in a Newtonian fluid. 2023. hal-04101291v2

HAL Id: hal-04101291

https://hal.science/hal-04101291v2

Preprint submitted on 28 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical and computational framework for moving and colliding rigid bodies in a Newtonian fluid

Céline Van Landeghem, Luca Berti, Vincent Chabannes, Agathe Chouippe, Laëtitia Giraldi, Yannick Hoarau, Christophe Prud'Homme

► To cite this version:

Céline Van Landeghem, Luca Berti, Vincent Chabannes, Agathe Chouippe, Laëtitia Giraldi, et al. Mathematical and computational framework for moving and colliding rigid bodies in a Newtonian fluid. Annals of Mathematical Sciences and Applications, 2024, pp.59-89. 10.4310/AMSA.2024.v9.n1.a2. hal-04101291v2

HAL Id: hal-04101291

https://hal.science/hal-04101291v2

Submitted on 28 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical and computation framework for moving and colliding rigid bodies in a Newtonian fluid

Céline Van Landeghem, Luca Berti, Vincent Chabannes, Christophe Prud'homme

*Cemosis, IRMA UMR 7501, CNRS, Université de Strasbourg, France*

Laetitia Giraldi

CALISTO team, INRIA, Université Côte d'Azur, France

Agathe Chouippe, Yannick Hoarau

ICube Laboratory, UMR 7357, Université de Strasbourg, France

## Abstract

We studied numerically the dynamics of colliding rigid bodies in a Newtonian fluid. The finite element method is used to solve the fluid-body interaction and the fluid motion is described in the Arbitrary-Lagrangian-Eulerian framework. To model the interactions between bodies, we consider a repulsive collision-avoidance model, defined by R. Glowinski in [1]. The main emphasis in this work is the generalization of this collision model to multiple rigid bodies of arbitrary shape. Our model first uses a narrow-band fast marching method to detect the set of colliding bodies. Then, collision forces and torques are computed for these bodies via a general expression, which does not depend on their shape. Numerical experiments examining the performance of the narrow-band fast marching method and the parallel execution of the collision algorithm are discussed. We validate our model with literature results and show various applications of colliding bodies in two and three dimensions. In these applications, the bodies either move due to gravity, a flow, or can actuate themselves. Finally, we present a tool to create arbitrary shaped bodies in complex already discretized fluid domains, enabling conforming body-fluid interface and allowing to perform simulations of fluid-body interactions with collision treatment in these realistic environments. All simulations are conducted with the Feel++ open source library.

**Keywords:** Fluid-structure interaction, rigid body motion, collision simulation, Feel++, MSC 2010 65M60, 74F10, 76M10, 70E99

## Introduction

Fluid flows laden with particles are common in industrial and biological processes, such as fluidization, the cell transport in arteries, or the simulated motion of articulated micro-swimmers in the human body. Due to the particle volume and the confined environments, these processes are characterized by inter-particle interactions. The modeling of these interactions, based on collision detection algorithms, and the computation of lubrication and collision forces, is challenging for arbitrary shaped particles and add further com-
Mathematical and computation framework for moving and
colliding rigid bodies in a Newtonian fluid

Céline Van Landeghem, Luca Berti, Vincent Chabannes,
Christophe Prud'homme

*Cemosis, IRMA UMR 7501, CNRS, Université de Strasbourg, France*

Laetitia Giraldi

CALISTO team, INRIA, Université Côte d'Azur, France

Agathe Chouippe, Yannick Hoarau

ICube Laboratory, UMR 7357, Université de Strasbourg, France

Abstract

We studied numerically the dynamics of colliding rigid bodies in a Newtonian fluid. The finite element method is used to solve the fluid-body interaction and the fluid motion is described in the Arbitrary-Lagrangian-Eulerian framework. To model the interactions between bodies, we consider a repulsive collision-avoidance model, defined by R. Glowinski in [1]. The main emphasis in this work is the generalization of this collision model to multiple rigid bodies of arbitrary shape. Our model first uses a narrow-band fast marching method to detect the set of colliding bodies. Then, collision forces and torques are computed for these bodies via a general expression, which does not depend on their shape. Numerical experiments examining the performance of the narrow-band fast marching method and the parallel execution of the collision algorithm are discussed. We validate our model with literature results and show various applications of colliding bodies in two and three dimensions. In these applications, the bodies either move due to gravity, a flow, or can actuate themselves. Finally, we present a tool to create arbitrary shaped bodies in complex already discretized fluid domains, enabling conforming body-fluid interface and allowing to perform simulations of fluid-body interactions with collision treatment in these realistic environments. All simulations are conducted with the Feel++ open source library.

**Keywords:** Fluid-structure interaction, rigid body motion, collision simulation, Feel++, MSC 2010 65M60, 74F10, 76M10, 70E99

**Introduction**

Fluid flows laden with particles are common in industrial and biological processes, such as fluidization, the cell transport in arteries, or the simulated motion of articulated micro-swimmers in the human body. Due to the particle volume and the confined environments, these processes are characterized by inter-particle interactions. The modeling of these interactions, based on collision detection algorithms, and the computation of lubrication and collision forces, is challenging for arbitrary shaped particles and add further com-
plexity to the coupled fluid-solid interaction problem.

In recent years, various approaches simulating solid-solid or solid-wall interactions have been developed and proposed in the literature. These approaches are based on collision and lubrication forces, added to the system when solving the fluid-solid inter-actions [2]. When the distance between solid surfaces is small, approximately the mesh size of the computational grid, the fluid placed between the solids is squeezed out and the hydrodynamics are not resolved. The resulting underestimation of the hydrodynamic forces is compensated by a short-range repulsive correction given analytically and called lubrication force. In [3], this lubrication force is defined by a constant expression and in [4] its magnitude is inversely proportional to the distance between the solid surfaces. In some approaches, the lubrication model allows avoiding the direct contact between the solid surfaces [1]. Such contact avoidance schemes are frequently used in literature to simulate the interactions between spherical bodies. In this paper, we will present the generalization of these schemes for complex shaped and articulated bodies. The contact avoidance models often require the adjustment of stiffness parameters, for which the opti-mal values are not known. Over- or underestimation of these values can falsify the results. To overcome this issue, the authors in [5] solve a minimization problem to determine the minimal force magnitude ensuring a fixed separation distance between the solids at each time step and in [6] the authors introduce the gluey particle method.

A collision force is applied, when the lubrication correction does not avoid the direct contact between the solids. Common collision forces are based on the soft-sphere or hard-sphere approach. When two solids are in direct contact or overlapping, the soft-sphere approach defines the collision force according to the geometrical characteristics of this contact [7]. The approach requires a small time step to resolve the collision process. In contrast, the hard-sphere approach [3] assumes that the collision force is impulsive. During the collision process, the velocities of the interacting solids change instantaneously at the time of contact. The magnitude of the new velocities depends on the pre-collision velocities and some physical parameters.

For large fluid systems, efficient collision detection algorithms are mandatory. These algorithms identify the pairs of solids that are actually interacting, before computing the collision forces only for the respective pairs, which reduces the computational costs. Most collision detection algorithms contain two phases, the broad-phase and the narrow-phase. First, the broad-phase uses spatial partitioning or sorting methods to determine the smallest possible set of neighbouring pairs that are likely to interact. In spatial partitioning, the fluid system is divided into regions and all solids assigned to the same region are considered as neighbours. The data structures employed for spatial partitioning and neighbour detection are lists [8], trees [9], adapted for large number of solids, or hash tables [10], efficient for solids with a wide range of sizes. The spatial sorting methods, as the Sort and Sweep algorithm [11], sort the solids in space to determine overlapping. Then, the narrow-phase detects whether these pairs are actually in contact by computing the distance between the solid surfaces. The distance computation being evident for spherical solids, it is a challenging and often expensive task for complex shaped or articulated bodies. In the literature, multiple analytical methods for specific geometries
plexity to the coupled fluid-solid interaction problem.

In recent years, various approaches simulating solid-solid or solid-wall interactions have been developed and proposed in the literature. These approaches are based on collision and lubrication forces, added to the system when solving the fluid-solid inter-actions [2]. When the distance between solid surfaces is small, approximately the mesh size of the computational grid, the fluid placed between the solids is squeezed out and the hydrodynamics are not resolved. The resulting underestimation of the hydrodynamic forces is compensated by a short-range repulsive correction given analytically and called lubrication force. In [3], this lubrication force is defined by a constant expression and in [4] its magnitude is inversely proportional to the distance between the solid surfaces. In some approaches, the lubrication model allows avoiding the direct contact between the solid surfaces [1]. Such contact avoidance schemes are frequently used in literature to simulate the interactions between spherical bodies. In this paper, we will present the generalization of these schemes for complex shaped and articulated bodies. The contact avoidance models often require the adjustment of stiffness parameters, for which the opti-mal values are not known. Over- or underestimation of these values can falsify the results. To overcome this issue, the authors in [5] solve a minimization problem to determine the minimal force magnitude ensuring a fixed separation distance between the solids at each time step and in [6] the authors introduce the gluey particle method.

A collision force is applied, when the lubrication correction does not avoid the direct contact between the solids. Common collision forces are based on the soft-sphere or hard-sphere approach. When two solids are in direct contact or overlapping, the soft-sphere approach defines the collision force according to the geometrical characteristics of this contact [7]. The approach requires a small time step to resolve the collision process. In contrast, the hard-sphere approach [3] assumes that the collision force is impulsive. During the collision process, the velocities of the interacting solids change instantaneously at the time of contact. The magnitude of the new velocities depends on the pre-collision velocities and some physical parameters.

For large fluid systems, efficient collision detection algorithms are mandatory. These algorithms identify the pairs of solids that are actually interacting, before computing the collision forces only for the respective pairs, which reduces the computational costs. Most collision detection algorithms contain two phases, the broad-phase and the narrow-phase. First, the broad-phase uses spatial partitioning or sorting methods to determine the smallest possible set of neighbouring pairs that are likely to interact. In spatial partitioning, the fluid system is divided into regions and all solids assigned to the same region are considered as neighbours. The data structures employed for spatial partitioning and neighbour detection are lists [8], trees [9], adapted for large number of solids, or hash tables [10], efficient for solids with a wide range of sizes. The spatial sorting methods, as the Sort and Sweep algorithm [11], sort the solids in space to determine overlapping. Then, the narrow-phase detects whether these pairs are actually in contact by computing the distance between the solid surfaces. The distance computation being evident for spherical solids, it is a challenging and often expensive task for complex shaped or articulated bodies. In the literature, multiple analytical methods for specific geometries
or general algorithms for arbitrary solid shapes are proposed [2]. Our model uses the fast marching method [12] to track the distance between any pair of arbitrary shaped bodies. Focusing on fluid systems containing up to one hundred solids, the implementation of a broad-phase algorithm is not necessary. However, to avoid the high computational costs associated to the fast marching method, we adapt a narrow-band approach, executable in parallel.

The paper is organized as follows. Section 1 details the coupled problem describing the fluid-solid interaction. The collision model including the collision detection algorithm based on the narrow-band fast marching method and the computation of the contact avoidance lubrication force is presented in section 2. The performance verification and the validation of the model are done in sections 3 and 4. Section 5 shows some applications of the collision model, including simulations of single complex shaped and articulated solids as well as spherical multi-bodies interactions. A summary and conclusions are given in section 6.

# 1 Rigid body moving in a fluid

This section briefly describes the coupled fluid-rigid body interaction problem. A detailed description is given in [13].

**The fluid model** The fluid model considered throughout the paper is incompressible and Newtonian, and the motion of the fluid domain is modelled using the Arbitrary-Lagrangian-Eulerian (ALE) formalism [14]. Let $\mathcal{F}_t \subseteq \mathbb{R}^d, d \in \{2, 3\}$, denote the domain occupied by the fluid at time $t$, where $t \in [0, T]$ and $T$ the final time of the simulation. Let $\mathcal{A}_t: \mathcal{F}_0 \to \mathcal{F}_t$ be the ALE map mapping the reference fluid domain $\mathcal{F}_0$ to the current domain $\mathcal{F}_t$, defined as $\mathcal{A}_t(X) = \mathcal{A}(t, X) = X + x(t, X)$ with $x(t, X)$ the displacement of the domain. Let $u: \mathcal{F}_t \times [0, T] \to \mathbb{R}^d$ and $p: \mathcal{F}_t \times [0, T] \to \mathbb{R}$ the fluid velocity and the hydrostatic pressure. Let $\rho_f$ and $\mu$ be the fluid density and dynamic viscosity, constant for the considered fluid model. Finally, let $\sigma \in \mathcal{M}_d(\mathbb{R})$ be the fluid stress tensor and $g \in \mathbb{R}^d$ the gravity acceleration. The ALE formulation of the Navier-Stokes equations partially decouples the geometric evolution of the fluid domain from that of the fluid continuum. Due to the change of frame, the ALE time derivative substitutes the Eulerian one $\partial_t u$ as $\partial_t u = \partial_t u|_A - (\partial_t x \cdot \nabla)u = \partial_t u|_x - (u_A \cdot \nabla)u$ in the momentum equation. The first term of the ALE derivative corresponds to the time variation of the fluid velocity as seen in the arbitrary frame, while the second contains the relative velocity between the fluid continuum and the new reference frame. The Navier-Stokes equations in the ALE frame are:

$$ \rho_f \partial_t u|_A + \rho_f \left( (u - u_A) \cdot \nabla \right) u = -\nabla \cdot \sigma + \rho_f g, \quad \text{in } \mathcal{F}_t, \tag{1} \\ \nabla \cdot u = 0, \quad \text{in } \mathcal{F}_t. $$

**The rigid body equations** A rigid body moving in a fluid is described by the motion of its center of mass, given by the Newton equation, and the rotation matrix between its local frame and the laboratory frame, defined by the Euler equation. In this paper, the
or general algorithms for arbitrary solid shapes are proposed [2]. Our model uses the fast marching method [12] to track the distance between any pair of arbitrary shaped bodies. Focusing on fluid systems containing up to one hundred solids, the implementation of a broad-phase algorithm is not necessary. However, to avoid the high computational costs associated to the fast marching method, we adapt a narrow-band approach, executable in parallel.

The paper is organized as follows. Section 1 details the coupled problem describing the fluid-solid interaction. The collision model including the collision detection algorithm based on the narrow-band fast marching method and the computation of the contact avoidance lubrication force is presented in section 2. The performance verification and the validation of the model are done in sections 3 and 4. Section 5 shows some applications of the collision model, including simulations of single complex shaped and articulated solids as well as spherical multi-bodies interactions. A summary and conclusions are given in section 6.

# 1 Rigid body moving in a fluid

This section briefly describes the coupled fluid-rigid body interaction problem. A detailed description is given in [13].

**The fluid model** The fluid model considered throughout the paper is incompressible and Newtonian, and the motion of the fluid domain is modelled using the Arbitrary-Lagrangian-Eulerian (ALE) formalism [14]. Let $\mathcal{F}_t \subseteq \mathbb{R}^d, d \in \{2, 3\}$, denote the domain occupied by the fluid at time $t$, where $t \in [0, T]$ and $T$ the final time of the simulation. Let $\mathcal{A}_t: \mathcal{F}_0 \to \mathcal{F}_t$ be the ALE map mapping the reference fluid domain $\mathcal{F}_0$ to the current domain $\mathcal{F}_t$, defined as $\mathcal{A}_t(X) = \mathcal{A}(t, X) = X + x(t, X)$ with $x(t, X)$ the displacement of the domain. Let $u: \mathcal{F}_t \times [0, T] \to \mathbb{R}^d$ and $p: \mathcal{F}_t \times [0, T] \to \mathbb{R}$ the fluid velocity and the hydrostatic pressure. Let $\rho_f$ and $\mu$ be the fluid density and dynamic viscosity, constant for the considered fluid model. Finally, let $\sigma \in \mathcal{M}_d(\mathbb{R})$ be the fluid stress tensor and $g \in \mathbb{R}^d$ the gravity acceleration. The ALE formulation of the Navier-Stokes equations partially decouples the geometric evolution of the fluid domain from that of the fluid continuum. Due to the change of frame, the ALE time derivative substitutes the Eulerian one $\partial_t u$ as $\partial_t u = \partial_t u|_A - (\partial_t x \cdot \nabla)u = \partial_t u|_x - (u_A \cdot \nabla)u$ in the momentum equation. The first term of the ALE derivative corresponds to the time variation of the fluid velocity as seen in the arbitrary frame, while the second contains the relative velocity between the fluid continuum and the new reference frame. The Navier-Stokes equations in the ALE frame are:

$$ \rho_f \partial_t u|_A + \rho_f \left( (u - u_A) \cdot \nabla \right) u = -\nabla \cdot \sigma + \rho_f g, \quad \text{in } \mathcal{F}_t, \tag{1} \\ \nabla \cdot u = 0, \quad \text{in } \mathcal{F}_t. $$

**The rigid body equations** A rigid body moving in a fluid is described by the motion of its center of mass, given by the Newton equation, and the rotation matrix between its local frame and the laboratory frame, defined by the Euler equation. In this paper, the
motion of a body is caused by fluid stresses and gravity, contact and lubrication forces that act as external forces $F_e$ and torques $T_e$. Let $S \subset \mathbb{R}^d$ be the domain occupied by the body, $\rho_S \in \mathbb{R}_{>0}$ its density and $m = \int_S \rho_S$ its mass. Let $U : [0, T] \to \mathbb{R}^d$ and $\omega : [0, T] \to \mathbb{R}^{d*}$, where $d^* = 1$ if $d=2$ or $d^* = 3$ if $d=3$, be the linear and angular velocity of the rigid body as seen from the laboratory frame. Let $x^{CM} = m^{-1} \int_S \rho_S x$ be its center of mass, $x^{CM} \in \mathbb{R}^d$, and $I = \int_S \rho_S(x - x^{CM}) \otimes (x - x^{CM})$ its inertia tensor, $I \in S_{++}^{d*}$ positive definite and symmetric. We chose to describe the rotation matrix $R(\theta) : \Theta \to SO(d)$ using Euler angles $\theta \in \Theta$, where $\Theta = [-\pi, \pi]$ if $d=2$, or $\Theta = [-\pi, \pi] \times [0, \pi] \times [0, \pi/2]$ if $d=3$. The Newton and Euler equations, describing the dynamics of a three-dimensional rigid body, are:

$$ m \frac{d}{dt} U = F_e - \int_{\partial S} [-pn + \mu(\nabla u + \nabla u^T)n], \\ \frac{d}{dt}(RIR^T \omega) = T_e - \int_{\partial S} [-pn + \mu(\nabla u + \nabla u^T)n] \times (x - x^{CM}), \quad (2) $$

where $n$ is the unit outward normal to $\partial S$ and one has:

$$ \frac{d}{dt} \theta_i = \omega_i, \quad \text{for } i \in \{x, y, z\}, \\ R = R_z(\theta_z) R_y(\theta_y) R_x(\theta_x), $$

where $R(\theta_i)$ denotes the rotation matrix around axis $i \in \{x, y, z\}$ of angle $\theta_i$. If $d=2$, $R(\theta)$ has the form:

$$ R(\theta) = \begin{bmatrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix}, $$

while in three dimensions:

$$ R_z(\theta_z) = \begin{bmatrix} \cos(\theta_z) & -\sin(\theta_z) & 0 \\ \sin(\theta_z) & \cos(\theta_z) & 0 \\ 0 & 0 & 1 \end{bmatrix}, R_y(\theta_y) = \begin{bmatrix} \cos(\theta_y) & 0 & \sin(\theta_y) \\ 0 & 1 & 0 \\ -\sin(\theta_y) & 0 & \cos(\theta_y) \end{bmatrix}, \\ R_x(\theta_x) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta_x) & -\sin(\theta_x) \\ 0 & \sin(\theta_x) & \cos(\theta_x) \end{bmatrix}. $$

**Fluid-solid interaction** The interaction between the rigid body and the fluid is ensured by the balance of stresses by the Newton and Euler equations (2) and a coupling condition at the interface, imposing the continuity of the velocities:

$$ u = U + \omega \times (x - x^{CM}) \quad \text{on } \partial \mathcal{F}_t \cap \partial S. \quad (3) $$

Finally, the definition of the $A_t(X)$ requires the computation of the displacement $x(t, X)$ in the reference domain. In order to find $x(t, X)$, we solve:

$$ \nabla \cdot ([1 + \tau(X)]\nabla_X x(t, X)) = 0 \quad \text{in } \mathcal{F}_0, \\ x(t, X) = g(t, X) \quad \text{in } \partial \mathcal{F}_0, \quad (4) $$

where $g(t, X) = \int_0^t U + \omega \times (X - X^{CM})$ is the rigid displacement of $S$ and $\tau(X)$ is a space-dependent coefficient, related to the volume of the simplexes in the triangulation.
motion of a body is caused by fluid stresses and gravity, contact and lubrication forces that act as external forces $F_e$ and torques $T_e$. Let $\mathcal{S} \subset \mathbb{R}^d$ be the domain occupied by the body, $\rho_S \in \mathbb{R}_{>0}$ its density and $m = \int_S \rho_S$ its mass. Let $U : [0, T] \to \mathbb{R}^d$ and $\omega : [0, T] \to \mathbb{R}^{d*}$, where $d^* = 1$ if $d=2$ or $d^* = 3$ if $d=3$, be the linear and angular velocity of the rigid body as seen from the laboratory frame. Let $x^{CM} = m^{-1} \int_S \rho_S x$ be its center of mass, $x^{CM} \in \mathbb{R}^d$, and $I = \int_S \rho_S(x - x^{CM}) \otimes (x - x^{CM})$ its inertia tensor, $I \in S_{++}^{d*}$ positive definite and symmetric. We chose to describe the rotation matrix $R(\theta) : \Theta \to SO(d)$ using Euler angles $\theta \in \Theta$, where $\Theta = [-\pi, \pi]$ if $d=2$, or $\Theta = [-\pi, \pi] \times [0, \pi] \times [0, \pi/2]$ if $d=3$. The Newton and Euler equations, describing the dynamics of a three-dimensional rigid body, are:

$$m \frac{d}{dt} U = F_e - \int_{\partial \mathcal{S}} -pn + \mu(\nabla u + \nabla u^T)n, \\ \frac{d}{dt}(RIR^T \omega) = T_e - \int_{\partial \mathcal{S}} [-pn + \mu(\nabla u + \nabla u^T)n] \times (x - x^{CM}), \quad (2)$$

where $n$ is the unit outward normal to $\partial \mathcal{S}$ and one has:

$$\frac{d}{dt} \theta_i = \omega_i, \quad \text{for } i \in \{x, y, z\}, \\ R = R_z(\theta_z) R_y(\theta_y) R_x(\theta_x),$$

where $R(\theta_i)$ denotes the rotation matrix around axis $i \in \{x, y, z\}$ of angle $\theta_i$. If $d=2$, $R(\theta)$ has the form:

$$R(\theta) = \begin{bmatrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix},$$

while in three dimensions:

$$R_z(\theta_z) = \begin{bmatrix} \cos(\theta_z) & -\sin(\theta_z) & 0 \\ \sin(\theta_z) & \cos(\theta_z) & 0 \\ 0 & 0 & 1 \end{bmatrix}, R_y(\theta_y) = \begin{bmatrix} \cos(\theta_y) & 0 & \sin(\theta_y) \\ 0 & 1 & 0 \\ -\sin(\theta_y) & 0 & \cos(\theta_y) \end{bmatrix}, \\ R_x(\theta_x) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta_x) & -\sin(\theta_x) \\ 0 & \sin(\theta_x) & \cos(\theta_x) \end{bmatrix}.$$

**Fluid-solid interaction** The interaction between the rigid body and the fluid is ensured by the balance of stresses by the Newton and Euler equations (2) and a coupling condition at the interface, imposing the continuity of the velocities:

$$u = U + \omega \times (x - x^{CM}) \quad \text{on } \partial \mathcal{F}_t \cap \partial \mathcal{S}. \quad (3)$$

Finally, the definition of the $\mathcal{A}_t(X)$ requires the computation of the displacement $x(t, X)$ in the reference domain. In order to find $x(t, X)$, we solve:

$$\nabla \cdot ([1 + \tau(X)]\nabla_X x(t, X)) = 0 \quad \text{in } \mathcal{F}_0, \\ x(t, X) = g(t, X) \quad \text{in } \partial \mathcal{F}_0, \quad (4)$$

where $g(t, X) = \int_0^t U + \omega \times (X - X^{CM})$ is the rigid displacement of $\mathcal{S}$ and $\tau(X)$ is a space-dependent coefficient, related to the volume of the simplexes in the triangulation.
**Numerical solution** The problem is discretized over a triangulation of the fluid domain. The discrete fluid velocity $u_h$ and pressure $p_h$ belong to the inf-sup stable Taylor-Hood space, where $P_2$ continuous finite elements are chosen for the velocity and $P_1$ continuous finite elements for the pressure. The discrete ALE map $A_h^t$ is also discretized using $P_1$ continuous finite elements. The coupling condition of the fluid-solid interaction is encoded in the finite element test spaces, as proposed in [15], which leads to the construction of an operator $P$ that couples the velocities at the interface, at the discrete level. In our computations, the interface is conforming and the solution of equation (4) handles small mesh deformations, with $\tau(X)$ a piecewise constant coefficient, defined on each element $e$ of the mesh as $\tau|_e = (1 - V_{min}/V_{max})/(V_e/V_{max})$, where $V_{max}$, $V_{min}$ and $V_e$ are the volumes of the largest, smallest and current element of the domain discretization [16]. In case of larger domain deformations, we perform remeshing and preserve the discretization of the interface at the same time.

## 2 Collision model

In this section, we will detail our collision model characterized by two phases: a collision detection algorithm and a lubrication model. We do not add a scheme defining direct contact between body surfaces, since our lubrication model avoids this kind of interaction. We will first describe the three types of rigid bodies for which we have developed our collision model. Then, we will detail the collision detection algorithm before defining the lubrication force.

### 2.1 Rigid body types

We adjust the collision model for three different types of rigid bodies. First, we consider spherical bodies in two and three dimensions. As collision models for spherical bodies are common in literature, we use this case for validation. Then, we adapt this first model for bodies of complex shape. In this paper, ellipses are used for the two-dimensional simulations and ellipsoids for simulations in three dimensions. At last, we consider articulated bodies, and in particular the three-sphere swimmer [17]. This type of swimmer is composed of three spheres of the same size that are connected by rods. The swimmer extends and retracts these rods to move. First the left sphere is retracted, then the right sphere. Finally, the left sphere is extended before extending the right sphere. This sequence of four movements result in a straight motion.

### 2.2 Collision detection

The first part of our model is a collision detection algorithm. Collision detection is an important task since it allows identifying the pairs of bodies that are actually interacting. Two bodies are interacting when the distance between their surfaces is smaller than the width of the collision zone, denoted $\rho$. In absence of detection phase, collision forces are computed for each pair of bodies present in the domain, independently of whether there will be an interaction. Most collision detection algorithms include two phases. First, the broad-phase identifies pairs of bodies in close neighborhood. Then, the distance between
**Numerical solution** The problem is discretized over a triangulation of the fluid domain. The discrete fluid velocity $u_h$ and pressure $p_h$ belong to the inf-sup stable Taylor-Hood space, where $P_2$ continuous finite elements are chosen for the velocity and $P_1$ continuous finite elements for the pressure. The discrete ALE map $A_h^t$ is also discretized using $P_1$ continuous finite elements. The coupling condition of the fluid-solid interaction is encoded in the finite element test spaces, as proposed in [15], which leads to the construction of an operator $P$ that couples the velocities at the interface, at the discrete level. In our computations, the interface is conforming and the solution of equation (4) handles small mesh deformations, with $\tau(X)$ a piecewise constant coefficient, defined on each element $e$ of the mesh as $\tau|_e = (1 - V_{min}/V_{max})/(V_e/V_{max})$, where $V_{max}$, $V_{min}$ and $V_e$ are the volumes of the largest, smallest and current element of the domain discretization [16]. In case of larger domain deformations, we perform remeshing and preserve the discretization of the interface at the same time.

## 2 Collision model

In this section, we will detail our collision model characterized by two phases: a collision detection algorithm and a lubrication model. We do not add a scheme defining direct contact between body surfaces, since our lubrication model avoids this kind of interaction. We will first describe the three types of rigid bodies for which we have developed our collision model. Then, we will detail the collision detection algorithm before defining the lubrication force.

### 2.1 Rigid body types

We adjust the collision model for three different types of rigid bodies. First, we consider spherical bodies in two and three dimensions. As collision models for spherical bodies are common in literature, we use this case for validation. Then, we adapt this first model for bodies of complex shape. In this paper, ellipses are used for the two-dimensional simulations and ellipsoids for simulations in three dimensions. At last, we consider articulated bodies, and in particular the three-sphere swimmer [17]. This type of swimmer is composed of three spheres of the same size that are connected by rods. The swimmer extends and retracts these rods to move. First the left sphere is retracted, then the right sphere. Finally, the left sphere is extended before extending the right sphere. This sequence of four movements result in a straight motion.

### 2.2 Collision detection

The first part of our model is a collision detection algorithm. Collision detection is an important task since it allows identifying the pairs of bodies that are actually interacting. Two bodies are interacting when the distance between their surfaces is smaller than the width of the collision zone, denoted $\rho$. In absence of detection phase, collision forces are computed for each pair of bodies present in the domain, independently of whether there will be an interaction. Most collision detection algorithms include two phases. First, the broad-phase identifies pairs of bodies in close neighborhood. Then, the distance between
the identified pairs is computed explicitly during the narrow-phase. According to this distance, one concludes if collision forces are applied. Our collision detection algorithm consists only of the narrow-phase. Thus, the distance between each pair of bodies is computed. Depending on the body shape, this is done in two different ways. We will first describe it for spherical and then for complex shaped or articulated bodies. Both methods require some parameters defined in a pre-processing step, which will also be described.

**Algorithm 1 Pre-process phase**

**Input:** Computational domain: mesh M

**Output:** Set of bodies identifiers: bodyIds, set of centers of mass: massCenters, set of bodies markers: bodyMarkers, set of radii: radii, set of centers of mass of imaginary bodies: ImagmassCenters, and the fluid marker: fluidMarker

**for** S in M.bodies **do**

bodyIds.append(S.id())

massCenters.append(S.massCenter())

bodyMarkers.append(S.boundary().name())

**if** type = spherical **then**

radii.append(S-radius())

ImagmassCenters.append(S.ImagmassCenters())

**end if**

**end for**

fluidMarker ← F_t.boundary().name()

**Spherical bodies** The distance $d_{ij}$ between the surfaces of two spherical bodies $S_i$ and $S_j$, $i \neq j$, is determined by first computing the distance between their mass centers $x_i^{CM}$ and $x_j^{CM}$ and then subtracting their radii $r_i$ and $r_j$:

$$d_{ij} = \|x_j^{CM} - x_i^{CM}\|_2 - r_i - r_j.$$

When this distance is smaller than a given parameter $\rho$, representing the width of the collision zone (detailed in subsection 2.3), then the identifiers $i, j$ and the distance $d_{ij}$ are stored in a collision map. Stored data is used during the phase where lubrication forces are computed.

The same formula can be used to get the distance $d_i$ between the surface of one body $S_i$ and the boundary of the fluid domain $\partial F_t$:

$$d_i = \|x_{i'}^{CM} - x_i^{CM}\|_2 - 2r_i.$$

where $x_{i'}^{CM}$ is the center of mass of the nearest imaginary body placed on the outside of the fluid boundary [18]. This imaginary body has same radius as the real body. Again, if $d_i$ is smaller than the width $\rho$, then the data $i$ and $d_i$ are stored. The parameters used in these two equations, i.e. body identifiers, mass centers and radii, are defined in a pre-processing step given by algorithm 1.

**Complex shaped or articulated bodies** When considering complex shaped or articulated bodies, then no explicit expression is provided to compute the distance between
the identified pairs is computed explicitly during the narrow-phase. According to this distance, one concludes if collision forces are applied. Our collision detection algorithm consists only of the narrow-phase. Thus, the distance between each pair of bodies is computed. Depending on the body shape, this is done in two different ways. We will first describe it for spherical and then for complex shaped or articulated bodies. Both methods require some parameters defined in a pre-processing step, which will also be described.

**Algorithm 1 Pre-process phase**

**Input:** Computational domain: mesh M

**Output:** Set of bodies identifiers: bodyIds, set of centers of mass: massCenters, set of bodies markers: bodyMarkers, set of radii: radii, set of centers of mass of imaginary bodies: ImagmassCenters, and the fluid marker: fluidMarker

**for** *S* in M.bodies **do**

bodyIds.append(*S*.id())

massCenters.append(*S*.massCenter())

bodyMarkers.append(*S*.boundary().name())

**if** type = spherical **then**

radii.append(*S*.radius())

ImagmassCenters.append(*S*.ImagmassCenters())

**end if**

**end for**

fluidMarker ← *F*<sub>*t*</sub>.boundary().name()

**Spherical bodies** The distance $d_{ij}$ between the surfaces of two spherical bodies $S_i$ and $S_j$, $i \neq j$, is determined by first computing the distance between their mass centers $x_i^{CM}$ and $x_j^{CM}$ and then subtracting their radii $r_i$ and $r_j$:

$$d_{ij} = ||x_j^{CM} - x_i^{CM}||_2 - r_i - r_j.$$

When this distance is smaller than a given parameter $\rho$, representing the width of the collision zone (detailed in subsection 2.3), then the identifiers $i, j$ and the distance $d_{ij}$ are stored in a collision map. Stored data is used during the phase where lubrication forces are computed.

The same formula can be used to get the distance $d_i$ between the surface of one body $S_i$ and the boundary of the fluid domain $\partial F_t$:

$$d_i = ||x_{i'}^{CM} - x_i^{CM}||_2 - 2r_i.$$

where $x_{i'}^{CM}$ is the center of mass of the nearest imaginary body placed on the outside of the fluid boundary [18]. This imaginary body has same radius as the real body. Again, if $d_i$ is smaller than the width $\rho$, then the data $i$ and $d_i$ are stored. The parameters used in these two equations, i.e. body identifiers, mass centers and radii, are defined in a pre-processing step given by algorithm 1.

**Complex shaped or articulated bodies** When considering complex shaped or articulated bodies, then no explicit expression is provided to compute the distance between
<caption>**Algorithm 2 Collision detection algorithm for spherical bodies**</caption>

**Input:** Output of pre-process phase
**Output:** Collision map
**for** i in bodyIds **do**
    **for** j in bodyIds, j > i **do**
        $d_{ij} \leftarrow \|\text{massCenters}_i - \text{massCenters}_j\|_2 - \text{radii}_i - \text{radii}_j$
        **if** $d_{ij} \le \rho$ **then**
            Store i, j and $d_{ij}$ in collision map
        **end if**
    **end for**
    **for** ImagCoord in ImagmassCenters_i **do**
        $d_i \leftarrow \|\text{massCenters}_i - \text{ImagCoord}\|_2 - 2\text{radii}_i$
        **if** $d_i \le \rho$ **then**
            $j \leftarrow -1$
            Store i, j, ImagCoord and $d_i$ in collision map
        **end if**
    **end for**
    **end for**

the surfaces. To deal with this issue, we use the fast marching method, based on a level set algorithm introduced by J.A Sethian [12]. By applying it to a body $S_i$, it returns the distance field $D_i$ from that body to the rest of the domain. Thus, the distance field $D_i$ is equal to zero at the body's boundary and has a positive value elsewhere. Given the distance fields $D_i$ and $D_j$, then the minimum distance between the surfaces $\partial S_i$ and $\partial S_j$ of the corresponding bodies is found by:

$$d_{ij} = \| \arg \min_{x \in \partial S_j} D_i(x) - \arg \min_{x \in \partial S_i} D_j(x) \|_2.$$

The minima $X_i = \arg \min_{x \in \partial S_i} D_j(x)$ and $X_j = \arg \min_{x \in \partial S_j} D_i(x)$ represent the contact points of the bodies $S_i$ and $S_j$, i.e. the coordinates of the surface points where the bodies will interact. To get the minimal distance between a body surface and the domain boundary, we use:

$$d_i = \| \arg \min_{x \in \partial \mathcal{F}_t} D_i(x) - \arg \min_{x \in \partial S_i} D_{\mathcal{F}_t}(x) \|_2,$$

where $D_{\mathcal{F}_t}$ the distance field obtained by applying the fast marching method to the domain boundary $\partial \mathcal{F}_t \setminus (\cup_i \partial S_i)$. The parameters that will be stored if $d_{ij} \le \rho$ or $d_i \le \rho$, are the identifiers, the minimal distance and the contact points.

When the fluid domain is large, the fast marching algorithm gets computationally expensive, especially in three dimensions. To accelerate the method, we develop the narrow-band approach, which computes the distance field only in a near neighborhood of the body boundary. This neighborhood is set to a predefined threshold $d_{max}$. As collision forces are only applied on a collision zone of width $\rho$, the threshold $d_{max}$ can be defined as function of $\rho$. Once the threshold distance $d_{max}$ is reached, the narrow-band approach assigns a default value $\delta$ to the distance field, corresponding to the maximum value reached:
<caption>**Algorithm 2** Collision detection algorithm for spherical bodies</caption>

**Input:** Output of pre-process phase
**Output:** Collision map
**for** i in bodyIds **do**
    **for** j in bodyIds, j > i **do**
        $d_{ij} \leftarrow \|\text{massCenters}_i - \text{massCenters}_j\|_2 - \text{radii}_i - \text{radii}_j$
        **if** $d_{ij} \le \rho$ **then**
            Store i, j and $d_{ij}$ in collision map
        **end if**
    **end for**
    **for** ImagCoord in ImagmassCenters_i **do**
        $d_i \leftarrow \|\text{massCenters}_i - \text{ImagCoord}\|_2 - 2\text{radii}_i$
        **if** $d_i \le \rho$ **then**
            $j \leftarrow -1$
            Store i, j, ImagCoord and $d_i$ in collision map
        **end if**
    **end for**
**end for**

the surfaces. To deal with this issue, we use the fast marching method, based on a level set algorithm introduced by J.A Sethian [12]. By applying it to a body $S_i$, it returns the distance field $D_i$ from that body to the rest of the domain. Thus, the distance field $D_i$ is equal to zero at the body's boundary and has a positive value elsewhere. Given the distance fields $D_i$ and $D_j$, then the minimum distance between the surfaces $\partial S_i$ and $\partial S_j$ of the corresponding bodies is found by:

$$d_{ij} = \| \arg \min_{x \in \partial S_j} D_i(x) - \arg \min_{x \in \partial S_i} D_j(x) \|_2.$$

The minima $X_i = \arg \min_{x \in \partial S_i} D_j(x)$ and $X_j = \arg \min_{x \in \partial S_j} D_i(x)$ represent the contact points of the bodies $S_i$ and $S_j$, i.e. the coordinates of the surface points where the bodies will interact. To get the minimal distance between a body surface and the domain boundary, we use:

$$d_i = \| \arg \min_{x \in \partial \mathcal{F}_t} D_i(x) - \arg \min_{x \in \partial S_i} D_{\mathcal{F}_t}(x) \|_2,$$

where $D_{\mathcal{F}_t}$ the distance field obtained by applying the fast marching method to the domain boundary $\partial \mathcal{F}_t \setminus (\cup_i \partial S_i)$. The parameters that will be stored if $d_{ij} \le \rho$ or $d_i \le \rho$, are the identifiers, the minimal distance and the contact points.

When the fluid domain is large, the fast marching algorithm gets computationally expensive, especially in three dimensions. To accelerate the method, we develop the narrow-band approach, which computes the distance field only in a near neighborhood of the body boundary. This neighborhood is set to a predefined threshold $d_{max}$. As collision forces are only applied on a collision zone of width $\rho$, the threshold $d_{max}$ can be defined as function of $\rho$. Once the threshold distance $d_{max}$ is reached, the narrow-band approach assigns a default value $\delta$ to the distance field, corresponding to the maximum value reached:
<caption>**Algorithm 3 Collision detection algorithm for complex shaped bodies**</caption>

**Input:** Output of pre-process phase
**Output:** Collision map
**for** *i* in bodyIds **do**
  *D*<sub>*i*</sub> ← narrow-band(bodyMarkers<sub>*i*</sub>, *d*<sub>max</sub> = 1.5*ρ*)
**end for**
*D*<sub>*F*<sub>*t*</sub></sub> ← narrow-band(fluidMarker, *d*<sub>max</sub> = 1.5*ρ*)
**for** *i* in bodyIds **do**
  **for** *j* in bodyIds, *j* > *i* **do**
    *d*<sub>*ij*</sub> ← ||*X*<sub>*j*</sub> − *X*<sub>*i*</sub>||<sub>2</sub>
    **if** *d*<sub>*ij*</sub> ≤ *ρ* **then**
      Store *i*, *j*, *d*<sub>*ij*</sub>, *X*<sub>*i*</sub> and *X*<sub>*j*</sub> in collision map
    **end if**
  **end for**
  *d*<sub>*i*</sub> ← ||*X*<sub>*F*<sub>*t*</sub></sub> − *X*<sub>*i*</sub>||<sub>2</sub>
  **if** *d*<sub>*i*</sub> ≤ *ρ* **then**
    *j* ← −1
    Store *i*, *j*, *d*<sub>*i*</sub>, *X*<sub>*i*</sub> and *X*<sub>*F*<sub>*t*</sub></sub> in collision map
  **end if**
**end for**

$$D_i^{NB} = \begin{cases} 0, & \text{on } \partial S_i, \\ D_i, & \text{for } D_i \le d_{max}, \\ \delta, & \text{elsewhere} \end{cases}$$

The narrow-band algorithm can be executed in parallel. A performance study of this approach is performed in section 3. In our implementation, the narrow-band fast marching function takes as parameter the boundary marker of the concerned body. Therefore, all boundary markers are stored during the pre-processing phase, presented in Algorithm 1.

When considering articulated bodies, the same formulas as for the complex shaped case are applied to each component of the body. Thus, when considering a pair of three-sphere swimmers, the distances between each sphere of one swimmer and all those of the second are computed. All steps of the collision detection algorithm for both methods are given by Algorithm 2 for spherical and by Algorithm 3 for complex shaped bodies. The computation of the imaginary bodies centers is only possible when the dimensions of the fluid domain, which must be rectangular, are known. For this reason, we have developed a more generic algorithm for spherical bodies, where the distance between a body and the fluid boundary is determined using the narrow-band approach. The initial algorithm can still be used for simulations of spherical bodies in rectangular domains, because its execution costs are lower.

## 2.3 Collision force

The second part of our collision model consists in the computation of collision forces. The lubrication model is based on a short-range repulsive force, a scheme introduced by
<caption>**Algorithm 3** Collision detection algorithm for complex shaped bodies</caption>

**Input:** Output of pre-process phase
**Output:** Collision map
**for** *i* in bodyIds **do**
  *D*<sub>*i*</sub> ← narrow-band(bodyMarkers<sub>*i*</sub>, *d*<sub>max</sub> = 1.5*ρ*)
**end for**
*D*<sub>*F*<sub>*t*</sub></sub> ← narrow-band(fluidMarker, *d*<sub>max</sub> = 1.5*ρ*)
**for** *i* in bodyIds **do**
  **for** *j* in bodyIds, *j* > *i* **do**
    *d*<sub>*ij*</sub> ← ||*X*<sub>*j*</sub> − *X*<sub>*i*</sub>||<sub>2</sub>
    **if** *d*<sub>*ij*</sub> ≤ *ρ* **then**
      Store *i*, *j*, *d*<sub>*ij*</sub>, *X*<sub>*i*</sub> and *X*<sub>*j*</sub> in collision map
    **end if**
  **end for**
  *d*<sub>*i*</sub> ← ||*X*<sub>*F*<sub>*t*</sub></sub> − *X*<sub>*i*</sub>||<sub>2</sub>
  **if** *d*<sub>*i*</sub> ≤ *ρ* **then**
    *j* ← −1
    Store *i*, *j*, *d*<sub>*i*</sub>, *X*<sub>*i*</sub> and *X*<sub>*F*<sub>*t*</sub></sub> in collision map
  **end if**
**end for**

$$D_i^{NB} = \begin{cases} 0, & \text{on } \partial S_i, \\ D_i, & \text{for } D_i \le d_{\text{max}}, \\ \delta, & \text{elsewhere} \end{cases}$$

The narrow-band algorithm can be executed in parallel. A performance study of this approach is performed in section 3. In our implementation, the narrow-band fast marching function takes as parameter the boundary marker of the concerned body. Therefore, all boundary markers are stored during the pre-processing phase, presented in Algorithm 1.

When considering articulated bodies, the same formulas as for the complex shaped case are applied to each component of the body. Thus, when considering a pair of three-sphere swimmers, the distances between each sphere of one swimmer and all those of the second are computed. All steps of the collision detection algorithm for both methods are given by Algorithm 2 for spherical and by Algorithm 3 for complex shaped bodies. The computation of the imaginary bodies centers is only possible when the dimensions of the fluid domain, which must be rectangular, are known. For this reason, we have developed a more generic algorithm for spherical bodies, where the distance between a body and the fluid boundary is determined using the narrow-band approach. The initial algorithm can still be used for simulations of spherical bodies in rectangular domains, because its execution costs are lower.

## 2.3 Collision force

The second part of our collision model consists in the computation of collision forces. The lubrication model is based on a short-range repulsive force, a scheme introduced by
R. Glowinski [1]. The force is activated when the distance between two body surfaces is less than the parameter $\rho$. This parameter represents the width of the collision zone, i.e. the zone where collision forces must be applied to prevent the bodies from overlapping and to avoid direct contact. The definition of the force is quite simple: it is parallel to the vector that connects the contact points of the bodies and its intensity increases as the distance decreases. We will first give the definition of this force for the simple case of spherical bodies. Then, we detail the required modifications, in particular the addition of its torque, for complex shaped and articulated bodies.

**Spherical bodies** For the repulsive force definition in the case of spherical bodies, we rely on articles [19, 20, 18]. The force to model the interaction between two bodies $S_i$ and $S_j$ is given by:

$$ \vec{F}_{ij} = \frac{1}{\epsilon} \overrightarrow{x_j^{CM} x_i^{CM}} (\rho - d_{ij})^2 \mathbf{1}_{d_{ij} \le \rho}. $$

In the same way, the interaction between a body $S_i$ and the domain boundary $\partial \mathcal{F}_t$ is defined by:

$$ \vec{F}_{i\mathcal{F}_t} = \frac{1}{\epsilon_{\mathcal{F}}} \overrightarrow{x_{i'}^{CM} x_i^{CM}} (\rho - d_i)^2 \mathbf{1}_{d_i \le \rho}. $$

Both equations contain a quadratic activation term. The vector connecting the mass centers gives the direction of the force, and the stiffness parameters $\epsilon$ and $\epsilon_{\mathcal{F}}$ describe its intensity for body-body and body-domain interactions. Finding the optimal values for these two parameters is not trivial, since their values depend on fluid and body properties. We refer to the indications given by article [18]: when the width of the collision zone is fixed to $\rho = 0.5h \sim 2.5h$, where $h$ is the mesh size, and the ratio between the body and fluid density is equal to 1, then one can suppose $\epsilon \approx h^2$ and $\epsilon_{\mathcal{F}} \approx \frac{\epsilon}{2}$. As already mentioned, these forces are only computed for pairs stored during collision detection phase, i.e. for pairs whose distance is less than $\rho$, in which cases the indicator functions $\mathbf{1}_{d_{ij} \le \rho} \equiv 1$ and $\mathbf{1}_{d_i \le \rho} \equiv 1$. The total repulsion force applied on one body $S_i$ is defined by:

$$ \vec{F}_i = \sum_{(i,j)|d_{ij} \le \rho} \vec{F}_{ij} + \sum_{i|d_i \le \rho} \vec{F}_{i\mathcal{F}_t}, $$

for all pairs $S_i - S_j$ and $S_i - \partial \mathcal{F}_t$ present in the collision map after detection algorithm. The total repulsion force $\vec{F}_i$ is added to the external forces $\vec{F}_e$ of the Newton equation in (2), describing the linear velocity of the body $S_i$.

**Complex shaped or articulated bodies** For the case of non-spherical bodies, the direction of the vector connecting the mass centers $x_i^{CM}$ and $x_j^{CM}$ no longer corresponds to that of the vector connecting the contact points $X_i$ and $X_j$. Thus, the equations are slightly modified:

$$ \vec{F}_{ij} = \frac{1}{\epsilon} \overrightarrow{X_j X_i} (\rho - d_{ij})^2 \mathbf{1}_{d_{ij} \le \rho}, \qquad \vec{F}_{i\mathcal{F}_t} = \frac{1}{\epsilon_{\mathcal{F}}} \overrightarrow{X_{i'} \mathcal{F}_t} (\rho - d_i)^2 \mathbf{1}_{d_i \le \rho}. $$

The values of the parameters $\epsilon$ and $\epsilon_{\mathcal{F}}$ have to be chosen smaller than those of the previous system. These modified expressions imply that the collision forces will lead to
R. Glowinski [1]. The force is activated when the distance between two body surfaces is less than the parameter $\rho$. This parameter represents the width of the collision zone, i.e. the zone where collision forces must be applied to prevent the bodies from overlapping and to avoid direct contact. The definition of the force is quite simple: it is parallel to the vector that connects the contact points of the bodies and its intensity increases as the distance decreases. We will first give the definition of this force for the simple case of spherical bodies. Then, we detail the required modifications, in particular the addition of its torque, for complex shaped and articulated bodies.

**Spherical bodies** For the repulsive force definition in the case of spherical bodies, we rely on articles [19, 20, 18]. The force to model the interaction between two bodies $S_i$ and $S_j$ is given by:

$$ \vec{F}_{ij} = \frac{1}{\epsilon} \overrightarrow{x_j^{CM} x_i^{CM}} (\rho - d_{ij})^2 \mathbf{1}_{d_{ij} \le \rho}. $$

In the same way, the interaction between a body $S_i$ and the domain boundary $\partial \mathcal{F}_t$ is defined by:

$$ \vec{F}_{i\mathcal{F}_t} = \frac{1}{\epsilon_{\mathcal{F}}} \overrightarrow{x_{i'}^{CM} x_i^{CM}} (\rho - d_i)^2 \mathbf{1}_{d_i \le \rho}. $$

Both equations contain a quadratic activation term. The vector connecting the mass centers gives the direction of the force, and the stiffness parameters $\epsilon$ and $\epsilon_{\mathcal{F}}$ describe its intensity for body-body and body-domain interactions. Finding the optimal values for these two parameters is not trivial, since their values depend on fluid and body properties. We refer to the indications given by article [18]: when the width of the collision zone is fixed to $\rho = 0.5h \sim 2.5h$, where $h$ is the mesh size, and the ratio between the body and fluid density is equal to 1, then one can suppose $\epsilon \approx h^2$ and $\epsilon_{\mathcal{F}} \approx \frac{\epsilon}{2}$. As already mentioned, these forces are only computed for pairs stored during collision detection phase, i.e. for pairs whose distance is less than $\rho$, in which cases the indicator functions $\mathbf{1}_{d_{ij} \le \rho} \equiv 1$ and $\mathbf{1}_{d_i \le \rho} \equiv 1$. The total repulsion force applied on one body $S_i$ is defined by:

$$ \vec{F}_i = \sum_{(i,j)|d_{ij} \le \rho} \vec{F}_{ij} + \sum_{i|d_i \le \rho} \vec{F}_{i\mathcal{F}_t}, $$

for all pairs $S_i - S_j$ and $S_i - \partial \mathcal{F}_t$ present in the collision map after detection algorithm. The total repulsion force $\vec{F}_i$ is added to the external forces $\vec{F}_e$ of the Newton equation in (2), describing the linear velocity of the body $S_i$.

**Complex shaped or articulated bodies** For the case of non-spherical bodies, the direction of the vector connecting the mass centers $x_i^{CM}$ and $x_j^{CM}$ no longer corresponds to that of the vector connecting the contact points $X_i$ and $X_j$. Thus, the equations are slightly modified:

$$ \vec{F}_{ij} = \frac{1}{\epsilon} \overrightarrow{X_j X_i} (\rho - d_{ij})^2 \mathbf{1}_{d_{ij} \le \rho}, \quad \vec{F}_{i\mathcal{F}_t} = \frac{1}{\epsilon_{\mathcal{F}}} \overrightarrow{X_{i'} \mathcal{F}_t} (\rho - d_i)^2 \mathbf{1}_{d_i \le \rho}. $$

The values of the parameters $\epsilon$ and $\epsilon_{\mathcal{F}}$ have to be chosen smaller than those of the previous system. These modified expressions imply that the collision forces will lead to
**Algorithm 4 Collision force algorithm for spherical bodies**

**Input:** Output of pre-process phase and collision map

**Output:** Set of collision forces: F, and set of collision torques: T

**for** pair in collision map **do**

if pair.j = -1 **then**

$$F_{i\mathcal{F}_t} \leftarrow \frac{1}{\epsilon} (\text{massCenters}_{\text{pair.i}} - \text{pair.ImagCoord})(\rho - \text{pair.d}_i)^2$$

$$F_{\text{pair.i}}.\text{add}(F_{i\mathcal{F}_t})$$

else

$$F_{ij} \leftarrow \frac{1}{\epsilon} (\text{massCenters}_{\text{pair.i}} - \text{massCenters}_{\text{pair.j}})(\rho - \text{pair.d}_{ij})^2$$

$$F_{\text{pair.i}}.\text{add}(F_{ij}), F_{\text{pair.j}}.\text{add}(-F_{ij})$$

end if

end for

body rotations. The torque to model these rotations is added to the external torques $T_e$ of the Euler equation in (2), and is defined by:

$$\vec{T}_i = -\overrightarrow{x_i^{CM}}\vec{X}_i \times \vec{F}_i.$$

This torque is computed for each body present in the collision map. The same approach is used for articulated bodies. The forces and torques are added for each element of the collision map. The algorithms given by Algorithm 4 and 5 illustrate the implementation of the lubrication model for spherical and complex shaped bodies.

**Algorithm 5 Collision force algorithm for complex shaped bodies**

**Input:** Output of pre-process phase and collision map

**Output:** Set of collision forces F, and set of collision torques T

**for** pair in collision map **do**

if pair.j = -1 **then**

$$F_{i\mathcal{F}_t} \leftarrow \frac{1}{\epsilon} (\text{pair}.X_i - \text{pair}.X_{\mathcal{F}_t})(\rho - \text{pair}.d_i)^2$$

$$F_{\text{pair.i}}.\text{add}(F_{i\mathcal{F}_t})$$

$$T_{\text{pair.i}}.\text{add}(-(\text{pair}.X_i - \text{massCenters}_{\text{pair.i}}).\text{cross}(F_{i\mathcal{F}_t}))$$

else

$$F_{ij} \leftarrow \frac{1}{\epsilon} (\text{pair}.X_i - \text{pair}.X_j)(\rho - \text{pair}.d_{ij})^2$$

$$F_{\text{pair.i}}.\text{add}(F_{ij}), F_{\text{pair.j}}.\text{add}(-F_{ij})$$

$$T_{\text{pair.i}}.\text{add}(-(\text{pair}.X_i - \text{massCenters}_{\text{pair.i}}).\text{cross}(F_{ij}))$$

$$T_{\text{pair.j}}.\text{add}(-(\text{pair}.X_j - \text{massCenters}_{\text{pair.j}}).\text{cross}(-F_{ij}))$$

end if

end for

# 3 Numerical experiments

## 3.1 Performance of narrow-band approach

The following test illustrates the performance of the narrow-band approach of the fast marching method, in two and three dimensions. The geometry consists in two spheres of same center but different radii. The radius of the inner and outer sphere is respectively set to $r_i = 0.1$ and $r_o = 2.0$. First, we use the fast marching method to determine the
**Algorithm 4 Collision force algorithm for spherical bodies**

**Input:** Output of pre-process phase and collision map

**Output:** Set of collision forces: F, and set of collision torques: T

**for pair in collision map do**

if pair.j = -1 **then**

$$F_{i\mathcal{F}_t} \leftarrow \frac{1}{\epsilon} (\text{massCenters}_{\text{pair.i}} - \text{pair.ImagCoord})(\rho - \text{pair.d}_i)^2$$

$$F_{\text{pair.i}}.\text{add}(F_{i\mathcal{F}_t})$$

else

$$F_{ij} \leftarrow \frac{1}{\epsilon} (\text{massCenters}_{\text{pair.i}} - \text{massCenters}_{\text{pair.j}})(\rho - \text{pair.d}_{ij})^2$$

$$F_{\text{pair.i}}.\text{add}(F_{ij}), F_{\text{pair.j}}.\text{add}(-F_{ij})$$

end if

end for

body rotations. The torque to model these rotations is added to the external torques $T_e$ of the Euler equation in (2), and is defined by:

$$\vec{T}_i = -\overrightarrow{x_i^{CM}}\vec{X}_i \times \vec{F}_i.$$

This torque is computed for each body present in the collision map. The same approach is used for articulated bodies. The forces and torques are added for each element of the collision map. The algorithms given by Algorithm 4 and 5 illustrate the implementation of the lubrication model for spherical and complex shaped bodies.

**Algorithm 5 Collision force algorithm for complex shaped bodies**

**Input:** Output of pre-process phase and collision map

**Output:** Set of collision forces F, and set of collision torques T

**for pair in collision map do**

if pair.j = -1 **then**

$$F_{i\mathcal{F}_t} \leftarrow \frac{1}{\epsilon} (\text{pair}.X_i - \text{pair}.X_{\mathcal{F}_t})(\rho - \text{pair}.d_i)^2$$

$$F_{\text{pair.i}}.\text{add}(F_{i\mathcal{F}_t})$$

$$T_{\text{pair.i}}.\text{add}(-(\text{pair}.X_i - \text{massCenters}_{\text{pair.i}}).\text{cross}(F_{i\mathcal{F}_t}))$$

else

$$F_{ij} \leftarrow \frac{1}{\epsilon} (\text{pair}.X_i - \text{pair}.X_j)(\rho - \text{pair}.d_{ij})^2$$

$$F_{\text{pair.i}}.\text{add}(F_{ij}), F_{\text{pair.j}}.\text{add}(-F_{ij})$$

$$T_{\text{pair.i}}.\text{add}(-(\text{pair}.X_i - \text{massCenters}_{\text{pair.i}}).\text{cross}(F_{ij}))$$

$$T_{\text{pair.j}}.\text{add}(-(\text{pair}.X_j - \text{massCenters}_{\text{pair.j}}).\text{cross}(-F_{ij}))$$

end if

end for

# 3 Numerical experiments

## 3.1 Performance of narrow-band approach

The following test illustrates the performance of the narrow-band approach of the fast marching method, in two and three dimensions. The geometry consists in two spheres of same center but different radii. The radius of the inner and outer sphere is respectively set to $r_i = 0.1$ and $r_o = 2.0$. First, we use the fast marching method to determine the
distance field from the inner sphere boundary to the outer sphere boundary. This is equivalent to perform the narrow-band approach on the inner sphere by setting the threshold distance $d_{max}$ to $d_{max} = r_o - r_i = 1.9$. Then, we continue to apply the narrow-band approach on the inner sphere but considering smaller thresholds. For the two-dimensional test the thresholds are set to $d_{max} = \{0.0625, 0.125, 0.25, 0.5, 1.0\}$, and we fix the mesh size to $h = 0.01$. The mesh contains $E_{2D} = 294918$ elements. In three dimensions, the mesh has $E_{3D} = 1385385$ elements, using a size equal to $h = 0.05$. The set of thresholds is $d_{max} = \{0.125, 0.25, 0.5, 1.0\}$.

The table 1 shows the execution time of the narrow-band approach for the different thresholds as well as the speed-up in two and three dimensions. We add the number of band elements $E_b$, i.e. the number of mesh elements in the zone where the distance field is actually computed, and the element ratio $\frac{E_{2D}}{E_b}$, where $n$ the dimension. It can be observed that the execution time decreases significantly when considering smaller thresholds. This test illustrates the importance of the narrow-band approach for our collision model.

<table><caption>Table 1: Results of performance test for the narrow-band approach in two and three dimensions. The table shows the execution time and speedup for different threshold values $d_{max}$.</caption><thead><tr><th>Threshold d<sub>max</sub></th><th>Band elements E<sub>b</sub></th><th>Element ratio E<sub>2D</sub>/E<sub>b</sub></th><th>Execution time in 2D</th><th>Speedup</th></tr></thead><tbody><tr><td>1.9</td><td>294918</td><td>1.00</td><td>1.959 s</td><td>/</td></tr><tr><td>1.0</td><td>89898</td><td>3.28</td><td>0.722 s</td><td>2.71</td></tr><tr><td>0.5</td><td>27082</td><td>10.89</td><td>0.411 s</td><td>4.76</td></tr><tr><td>0.25</td><td>9430</td><td>31.27</td><td>0.316 s</td><td>6.20</td></tr><tr><td>0.125</td><td>4104</td><td>71.86</td><td>0.276 s</td><td>7.09</td></tr><tr><td>0.0625</td><td>2262</td><td>130.38</td><td>0.261 s</td><td>7.50</td></tr></tbody><thead><tr><th>Threshold d<sub>max</sub></th><th>Band elements E<sub>b</sub></th><th>Element ratio E<sub>3D</sub>/E<sub>b</sub></th><th>Execution time in 3D</th><th>Speedup</th></tr></thead><tbody><tr><td>1.9</td><td>1385385</td><td>1.00</td><td>22.851 s</td><td>/</td></tr><tr><td>1.0</td><td>254365</td><td>5.45</td><td>5.869 s</td><td>3.89</td></tr><tr><td>0.5</td><td>43784</td><td>31.64</td><td>2.368 s</td><td>9.65</td></tr><tr><td>0.25</td><td>10238</td><td>135.32</td><td>1.815 s</td><td>12.59</td></tr><tr><td>0.125</td><td>4537</td><td>305.35</td><td>1.759 s</td><td>12.99</td></tr></tbody></table>

## 3.2 Performance of the parallel collision algorithm

The aim of this second numerical experiment is to observe the performance of the parallel implementation of the collision algorithm using the narrow-band fast marching method as well as to analyze the execution time proportion attributed to the collision model during the fluid-solid interaction resolution. The choice of the preconditioners, a direct LU solver in two dimensions and a block preconditioner in 3D, are detailed in [13]. We consider different numbers of spherical bodies distributed in a rectangular domain filled with a steady fluid. No gravity is applied and the magnitude of the collision forces is chosen close to zero, such that the bodies do not move during the simulation. Each
distance field from the inner sphere boundary to the outer sphere boundary. This is equivalent to perform the narrow-band approach on the inner sphere by setting the threshold distance $d_{max}$ to $d_{max} = r_o - r_i = 1.9$. Then, we continue to apply the narrow-band approach on the inner sphere but considering smaller thresholds. For the two-dimensional test the thresholds are set to $d_{max} = \{0.0625, 0.125, 0.25, 0.5, 1.0\}$, and we fix the mesh size to $h = 0.01$. The mesh contains $E_{2D} = 294918$ elements. In three dimensions, the mesh has $E_{3D} = 1385385$ elements, using a size equal to $h = 0.05$. The set of thresholds is $d_{max} = \{0.125, 0.25, 0.5, 1.0\}$.

The table 1 shows the execution time of the narrow-band approach for the different thresholds as well as the speed-up in two and three dimensions. We add the number of band elements $E_b$, i.e. the number of mesh elements in the zone where the distance field is actually computed, and the element ratio $\frac{E_{2D}}{E_b}$, where $n$ the dimension. It can be observed that the execution time decreases significantly when considering smaller thresholds. This test illustrates the importance of the narrow-band approach for our collision model.

<table><caption>Table 1: Results of performance test for the narrow-band approach in two and three dimensions. The table shows the execution time and speedup for different threshold values $d_{max}$.</caption><thead><tr><th>Threshold d<sub>max</sub></th><th>Band elements E<sub>b</sub></th><th>Element ratio <sup>E<sub>2D</sub></sup>&frasl;<sub>E<sub>b</sub></sub></th><th>Execution time in 2D</th><th>Speedup</th></tr></thead><tbody><tr><td>1.9</td><td>294918</td><td>1.00</td><td>1.959 s</td><td>/</td></tr><tr><td>1.0</td><td>89898</td><td>3.28</td><td>0.722 s</td><td>2.71</td></tr><tr><td>0.5</td><td>27082</td><td>10.89</td><td>0.411 s</td><td>4.76</td></tr><tr><td>0.25</td><td>9430</td><td>31.27</td><td>0.316 s</td><td>6.20</td></tr><tr><td>0.125</td><td>4104</td><td>71.86</td><td>0.276 s</td><td>7.09</td></tr><tr><td>0.0625</td><td>2262</td><td>130.38</td><td>0.261 s</td><td>7.50</td></tr></tbody><thead><tr><th>Threshold d<sub>max</sub></th><th>Band elements E<sub>b</sub></th><th>Element ratio <sup>E<sub>3D</sub></sup>&frasl;<sub>E<sub>b</sub></sub></th><th>Execution time in 3D</th><th>Speedup</th></tr></thead><tbody><tr><td>1.9</td><td>1385385</td><td>1.00</td><td>22.851 s</td><td>/</td></tr><tr><td>1.0</td><td>254365</td><td>5.45</td><td>5.869 s</td><td>3.89</td></tr><tr><td>0.5</td><td>43784</td><td>31.64</td><td>2.368 s</td><td>9.65</td></tr><tr><td>0.25</td><td>10238</td><td>135.32</td><td>1.815 s</td><td>12.59</td></tr><tr><td>0.125</td><td>4537</td><td>305.35</td><td>1.759 s</td><td>12.99</td></tr></tbody></table>

## 3.2 Performance of the parallel collision algorithm

The aim of this second numerical experiment is to observe the performance of the parallel implementation of the collision algorithm using the narrow-band fast marching method as well as to analyze the execution time proportion attributed to the collision model during the fluid-solid interaction resolution. The choice of the preconditioners, a direct LU solver in two dimensions and a block preconditioner in 3D, are detailed in [13]. We consider different numbers of spherical bodies distributed in a rectangular domain filled with a steady fluid. No gravity is applied and the magnitude of the collision forces is chosen close to zero, such that the bodies do not move during the simulation. Each
domain configuration is simulated for ten time iterations in sequential and in parallel, using different numbers of processors. For each case various information are displayed: the number of bodies, the number of mesh nodes, the number of body-body and body-wall interactions during each iteration as well as the speed-up of the collision algorithm (pre-process phase, collision detection algorithm, and computation of collision force) after the ten iterations. This speed-up $\frac{T_{np1}}{T_{npN}}$ is given by the ratio between the sequential execution time and the parallel execution time, when running the simulation on N processors. The results are given in tables 2 and 3.

<table><caption>Table 2: Results of the parallelization performance of the collision algorithm in two dimensions. The table shows the speed-up of the execution time when running the simulations in parallel.</caption><thead><tr><th>Bodies</th><th>Mesh nodes</th><th>Body-body</th><th>Body-wall</th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np4</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np16</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np24</sub></sub></th></tr></thead><tbody><tr><td>1</td><td>16462</td><td>0</td><td>1</td><td>3.14</td><td>8.22</td><td>11.04</td></tr><tr><td>25</td><td>17780</td><td>0</td><td>16</td><td>1.88</td><td>2.06</td><td>2.12</td></tr><tr><td>49</td><td>18284</td><td>28</td><td>24</td><td>1.90</td><td>2.08</td><td>2.15</td></tr><tr><td>81</td><td>18813</td><td>72</td><td>32</td><td>1.92</td><td>2.09</td><td>2.13</td></tr><tr><td>100</td><td>19139</td><td>100</td><td>36</td><td>1.81</td><td>2.01</td><td>2.10</td></tr></tbody></table>

<table><caption>Table 3: Parallelization performance in three dimensions. The sequential execution time is compared to the execution time of the simulation when running it in parallel using different numbers of processors.</caption><thead><tr><th>Bodies</th><th>Mesh nodes</th><th>Body-body</th><th>Body-wall</th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np4</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np16</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np24</sub></sub></th></tr></thead><tbody><tr><td>1</td><td>17073</td><td>0</td><td>1</td><td>2.23</td><td>6.49</td><td>7.99</td></tr><tr><td>27</td><td>18213</td><td>0</td><td>26</td><td>1.94</td><td>3.99</td><td>5.03</td></tr><tr><td>64</td><td>18789</td><td>0</td><td>56</td><td>1.98</td><td>3.98</td><td>5.22</td></tr><tr><td>125</td><td>20301</td><td>70</td><td>98</td><td>1.99</td><td>3.59</td><td>5.13</td></tr></tbody></table>

One can observe that the parallel execution reduces the execution times; the speed-up increases with an increasing number of processors. For multiple interactions, the speed-up remains almost constant independently of the number of bodies. The execution time spent in collision detection represents the majority of the total execution time of the algorithm. For the simulations conducted in this section, the execution time of the detection phase represents 99% of the total collision time. The difference in speed-ups between the two and three-dimensional case could be explained by the larger execution time taken by the distance computation, since the number of node-node connections in 3D is, on average, larger than in 2D. Furthermore, the execution time of the two-dimensional collision model represents up to 20% of the total time for the fluid-solid interaction resolution, when considering multiple interactions and using a direct solver. In three dimensions, using a block preconditioner, the proportion associated to the collision algorithm does not exceed 5%.
domain configuration is simulated for ten time iterations in sequential and in parallel, using different numbers of processors. For each case various information are displayed: the number of bodies, the number of mesh nodes, the number of body-body and body-wall interactions during each iteration as well as the speed-up of the collision algorithm (pre-process phase, collision detection algorithm, and computation of collision force) after the ten iterations. This speed-up $\frac{T_{np1}}{T_{npN}}$ is given by the ratio between the sequential execution time and the parallel execution time, when running the simulation on N processors. The results are given in tables 2 and 3.

<table><caption>Table 2: Results of the parallelization performance of the collision algorithm in two dimensions. The table shows the speed-up of the execution time when running the simulations in parallel.</caption><thead><tr><th>Bodies</th><th>Mesh nodes</th><th>Body-body</th><th>Body-wall</th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np4</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np16</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np24</sub></sub></th></tr></thead><tbody><tr><td>1</td><td>16462</td><td>0</td><td>1</td><td>3.14</td><td>8.22</td><td>11.04</td></tr><tr><td>25</td><td>17780</td><td>0</td><td>16</td><td>1.88</td><td>2.06</td><td>2.12</td></tr><tr><td>49</td><td>18284</td><td>28</td><td>24</td><td>1.90</td><td>2.08</td><td>2.15</td></tr><tr><td>81</td><td>18813</td><td>72</td><td>32</td><td>1.92</td><td>2.09</td><td>2.13</td></tr><tr><td>100</td><td>19139</td><td>100</td><td>36</td><td>1.81</td><td>2.01</td><td>2.10</td></tr></tbody></table>

<table><caption>Table 3: Parallelization performance in three dimensions. The sequential execution time is compared to the execution time of the simulation when running it in parallel using different numbers of processors.</caption><thead><tr><th>Bodies</th><th>Mesh nodes</th><th>Body-body</th><th>Body-wall</th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np4</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np16</sub></sub></th><th><sup>T<sub>np1</sub></sup>/<sub>T<sub>np24</sub></sub></th></tr></thead><tbody><tr><td>1</td><td>17073</td><td>0</td><td>1</td><td>2.23</td><td>6.49</td><td>7.99</td></tr><tr><td>27</td><td>18213</td><td>0</td><td>26</td><td>1.94</td><td>3.99</td><td>5.03</td></tr><tr><td>64</td><td>18789</td><td>0</td><td>56</td><td>1.98</td><td>3.98</td><td>5.22</td></tr><tr><td>125</td><td>20301</td><td>70</td><td>98</td><td>1.99</td><td>3.59</td><td>5.13</td></tr></tbody></table>

One can observe that the parallel execution reduces the execution times; the speed-up increases with an increasing number of processors. For multiple interactions, the speed-up remains almost constant independently of the number of bodies. The execution time spent in collision detection represents the majority of the total execution time of the algorithm. For the simulations conducted in this section, the execution time of the detection phase represents 99% of the total collision time. The difference in speed-ups between the two and three-dimensional case could be explained by the larger execution time taken by the distance computation, since the number of node-node connections in 3D is, on average, larger than in 2D. Furthermore, the execution time of the two-dimensional collision model represents up to 20% of the total time for the fluid-solid interaction resolution, when considering multiple interactions and using a direct solver. In three dimensions, using a block preconditioner, the proportion associated to the collision algorithm does not exceed 5%.
# 4 Validation

To validate our model we consider the motion of a circular solid in an incompressible Newtonian fluid in two and three dimensions.

For the validation in 2D, we consider a rigid body of radius $r = 0.125\text{cm}$ and density $\rho_s = 1.25\frac{\text{g}}{\text{cm}^3}$ is initially located at (1cm, 4cm) in a channel of width 2cm and height 6cm. Under the effect of gravity $g = 981\frac{\text{cm}}{\text{s}^2}$ the particle is falling towards the bottom of a channel filled with a fluid of density $\rho_f = 1\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 0.1\frac{\text{g}}{\text{cm s}}$. Both fluid and particle are at rest at time $t = 0.0\text{s}$. The mesh size and time step used in this simulation are respectively fixed at $h = 0.01\text{cm}$ and $\Delta t = 0.001\text{s}$. The width of the collision zone is set to $\rho = 0.015\text{cm}$. One can use a stiffness parameter for particle-fluid domain interaction $\epsilon_F < \frac{h^2}{2}$, for the results shown in the figures we used $\epsilon_F = 5 * 10^{-6}$. In figure 1, we plot the time evolution of four quantities and compare them to literature results Wan and Turek [18], Wang, Guo and Mi [19]. These quantities are the vertical coordinate of the center of mass $y^{cm}$, the vertical translational velocity $v_y$, the Reynolds number $Re$ and the translational kinetic energy $E_t$ defined by:

$$Re = \frac{2r\rho_s\sqrt{v_x^2 + v_y^2}}{\mu} \quad \text{and} \quad E_t = 0.5\pi r^2 \rho_s(v_x^2 + v_y^2),$$

where $v_x$ the horizontal translational velocity. One can observe that all results are in good agreement. The small differences before and after collision can be explained by different definitions of collision parameters and numerical methods.

We perform the same simulation for the 3D validation. Initially, a sphere of radius $r = 0.75\text{cm}$ and density $\rho_s = 1.12\frac{\text{g}}{\text{cm}^3}$ is located at (0cm, 12.75cm, 0cm) in a computational box of dimensions $[-4\text{cm}, 4\text{cm}] \times [0\text{cm}, 15\text{cm}] \times [-4\text{cm}, 4\text{cm}]$. The box is filled with a fluid of density $\rho_f = 0.962\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 1.13\frac{\text{g}}{\text{cm s}}$. The mesh size and the time step are set to $h = 0.1\text{cm}$ and $\Delta t = 0.01\text{s}$. Regarding the lubrication force parameters, we fixed the width of the collision zone to $\rho = 0.2625\text{cm}$ and the stiffness parameter to $\epsilon_F = 5 * 10^{-6}$. The comparison of the present results to literature results [21] is given in figure 2. We obtain the same results for the vertical coordinate of the center of mass and the vertical translational velocity.

# 5 Applications

## 5.1 Isolated object

In this subsection we simulate the interaction between a complex shaped body falling under the effect of gravity and the boundary of the computational domain filled with an incompressible Newtonian fluid in two and three dimensions.

For the first case, we consider a two-dimensional ellipse located in a channel of width $L = \frac{16}{130}\text{cm}$ and infinite length. The channel contains a fluid of density $\rho_f = 1.195\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 0.305\frac{\text{g}}{\text{cm s}}$. The ellipse's long and short axis are respectively set to
# 4 Validation

To validate our model we consider the motion of a circular solid in an incompressible Newtonian fluid in two and three dimensions.

For the validation in 2D, we consider a rigid body of radius $r = 0.125\text{cm}$ and density $\rho_s = 1.25\frac{\text{g}}{\text{cm}^3}$ is initially located at (1cm, 4cm) in a channel of width 2cm and height 6cm. Under the effect of gravity $g = 981\frac{\text{cm}}{\text{s}^2}$ the particle is falling towards the bottom of a channel filled with a fluid of density $\rho_f = 1\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 0.1\frac{\text{g}}{\text{cm s}}$. Both fluid and particle are at rest at time $t = 0.0\text{s}$. The mesh size and time step used in this simulation are respectively fixed at $h = 0.01\text{cm}$ and $\Delta t = 0.001\text{s}$. The width of the collision zone is set to $\rho = 0.015\text{cm}$. One can use a stiffness parameter for particle-fluid domain interaction $\epsilon_F < \frac{h^2}{2}$, for the results shown in the figures we used $\epsilon_F = 5 * 10^{-6}$. In figure 1, we plot the time evolution of four quantities and compare them to literature results Wan and Turek [18], Wang, Guo and Mi [19]. These quantities are the vertical coordinate of the center of mass $y^{cm}$, the vertical translational velocity $v_y$, the Reynolds number $Re$ and the translational kinetic energy $E_t$ defined by:

$$Re = \frac{2r\rho_s\sqrt{v_x^2 + v_y^2}}{\mu} \quad \text{and} \quad E_t = 0.5\pi r^2 \rho_s(v_x^2 + v_y^2),$$

where $v_x$ the horizontal translational velocity. One can observe that all results are in good agreement. The small differences before and after collision can be explained by different definitions of collision parameters and numerical methods.

We perform the same simulation for the 3D validation. Initially, a sphere of radius $r = 0.75\text{cm}$ and density $\rho_s = 1.12\frac{\text{g}}{\text{cm}^3}$ is located at (0cm, 12.75cm, 0cm) in a computational box of dimensions $[-4\text{cm}, 4\text{cm}] \times [0\text{cm}, 15\text{cm}] \times [-4\text{cm}, 4\text{cm}]$. The box is filled with a fluid of density $\rho_f = 0.962\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 1.13\frac{\text{g}}{\text{cm s}}$. The mesh size and the time step are set to $h = 0.1\text{cm}$ and $\Delta t = 0.01\text{s}$. Regarding the lubrication force parameters, we fixed the width of the collision zone to $\rho = 0.2625\text{cm}$ and the stiffness parameter to $\epsilon_F = 5 * 10^{-6}$. The comparison of the present results to literature results [21] is given in figure 2. We obtain the same results for the vertical coordinate of the center of mass and the vertical translational velocity.

# 5 Applications

## 5.1 Isolated object

In this subsection we simulate the interaction between a complex shaped body falling under the effect of gravity and the boundary of the computational domain filled with an incompressible Newtonian fluid in two and three dimensions.

For the first case, we consider a two-dimensional ellipse located in a channel of width $L = \frac{16}{130}\text{cm}$ and infinite length. The channel contains a fluid of density $\rho_f = 1.195\frac{\text{g}}{\text{cm}^3}$ and viscosity $\mu = 0.305\frac{\text{g}}{\text{cm s}}$. The ellipse's long and short axis are respectively set to
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the vertical position of mass center, vertical translational velocity, Reynolds number and translational kinetic energy. The red line are the present results, the green dotted and blue dotted lines correspond to results respectively taken from [19] and [18]. The streamlines and body position at four different time instants are represented.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the vertical position of mass center, vertical translational velocity, Reynolds number and translational kinetic energy. The red line are the present results, the green dotted and blue dotted lines correspond to results respectively taken from [19] and [18]. The streamlines and body position at four different time instants are represented.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison of the vertical position of the center of mass and vertical translational velocity observed in present results, red line, to literature results [21], dotted line.</figcaption></figure>

$a = 0.1\text{cm}$ and $b = 0.05\text{cm}$. Its density is equal to $\rho_S = 1.35\frac{g}{\text{cm}^3}$. Due to its initial orientation, set to $\theta = \frac{\pi}{3}$ with respect to the horizontal axis, the ellipse will collide with the right and left walls. To model these collisions the length of the lubrication zone is fixed at $\rho = 0.05\text{cm}$ and a stiffness parameter of $\epsilon_F = 10^{-5}$ is used. The simulation is performed using a mesh size $h = 0.001\text{cm}$ and a time step $\Delta t = 0.001\text{s}$. Figure 3a shows the trajectory of the ellipse comparing it to literature [22]. This trajectory allows to observe that the ellipse first performs oscillations between the right and left wall before reaching a horizontal position, at the instant given by the black point. Then, it performs rotations close to a single wall. In present results, these rotations are taking place near the left wall, but in the reference article, they happen near the right wall. Given that the ellipse is in a horizontal position at the beginning of these rotations, any small difference in the collision model or the numerical resolution techniques can explain this difference. To show that the present results remain close to literature, we consider in figure 3b the symmetric trajectory for the rotations of [22]. This test case allows us to validate our implementation for the case of arbitrarily shaped bodies.

The second case represents the simulation of a three-dimensional ellipsoid of density $\rho_S = 1.25\frac{g}{\text{cm}^3}$ and axis of length $a = 0.4\text{cm}$ and $b = c = 0.2\text{cm}$ falling in a rectangular computational domain filled with a fluid of density $\rho_f = 1\frac{g}{\text{cm}^3}$ and viscosity $\mu = 0.01\frac{g}{\text{cm s}}$. The dimensions of the domain are set to $[0\text{cm}, 1\text{cm}] \times [0\text{cm}, 8\text{cm}] \times [0\text{cm}, 0.4\text{cm}]$ and at $t = 0\text{s}$ the ellipsoid is located at $(0.5\text{cm}, 6\text{cm}, 0.2\text{cm})$ with its long axis oriented in vertical direction. The mesh size and time step of this simulation are $h = 0.0125\text{cm}$ and $\Delta t = 0.001\text{s}$. The interactions between the ellipsoid and the boundary of the computational domain are modelled by defining a collision zone of width $\rho = 0.04\text{cm}$ and a stiffness parameter $\epsilon_F = 3*10^{-6}$. The evolution of the horizontal position of the sphere is shown in figure 4, and is compared to literature results [23]. The trajectories before and
<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison of the vertical position of the center of mass and vertical translational velocity observed in present results, red line, to literature results [21], dotted line.</figcaption></figure>

$a = 0.1\text{cm}$ and $b = 0.05\text{cm}$. Its density is equal to $\rho_S = 1.35\frac{g}{\text{cm}^3}$. Due to its initial orientation, set to $\theta = \frac{\pi}{3}$ with respect to the horizontal axis, the ellipse will collide with the right and left walls. To model these collisions the length of the lubrication zone is fixed at $\rho = 0.05\text{cm}$ and a stiffness parameter of $\epsilon_F = 10^{-5}$ is used. The simulation is performed using a mesh size $h = 0.001\text{cm}$ and a time step $\Delta t = 0.001\text{s}$. Figure 3a shows the trajectory of the ellipse comparing it to literature [22]. This trajectory allows to observe that the ellipse first performs oscillations between the right and left wall before reaching a horizontal position, at the instant given by the black point. Then, it performs rotations close to a single wall. In present results, these rotations are taking place near the left wall, but in the reference article, they happen near the right wall. Given that the ellipse is in a horizontal position at the beginning of these rotations, any small difference in the collision model or the numerical resolution techniques can explain this difference. To show that the present results remain close to literature, we consider in figure 3b the symmetric trajectory for the rotations of [22]. This test case allows us to validate our implementation for the case of arbitrarily shaped bodies.

The second case represents the simulation of a three-dimensional ellipsoid of density $\rho_S = 1.25\frac{g}{\text{cm}^3}$ and axis of length $a = 0.4\text{cm}$ and $b = c = 0.2\text{cm}$ falling in a rectangular computational domain filled with a fluid of density $\rho_f = 1\frac{g}{\text{cm}^3}$ and viscosity $\mu = 0.01\frac{g}{\text{cm s}}$. The dimensions of the domain are set to $[0\text{cm}, 1\text{cm}] \times [0\text{cm}, 8\text{cm}] \times [0\text{cm}, 0.4\text{cm}]$ and at $t = 0\text{s}$ the ellipsoid is located at $(0.5\text{cm}, 6\text{cm}, 0.2\text{cm})$ with its long axis oriented in vertical direction. The mesh size and time step of this simulation are $h = 0.0125\text{cm}$ and $\Delta t = 0.001\text{s}$. The interactions between the ellipsoid and the boundary of the computational domain are modelled by defining a collision zone of width $\rho = 0.04\text{cm}$ and a stiffness parameter $\epsilon_F = 3*10^{-6}$. The evolution of the horizontal position of the sphere is shown in figure 4, and is compared to literature results [23]. The trajectories before and
<figure><img src="image_4.png" /><figcaption>Figure 3: The left graph plots the trajectory of the ellipse center of mass. The red line corresponds to the present results and the blue dotted line represents the literature results [22]. The right graph shows the same trajectory but compares it to rotations happening on the other wall.</figcaption></figure>

after the first interaction between the ellipsoid and the boundary of the computational
domain, at time $t \approx 0.45\text{s}$, are in good agreement.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the horizontal position of the center of the ellipsoid. The dotted line shows the literature results [23] and the red line represents the present results.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: The left graph plots the trajectory of the ellipse center of mass. The red line corresponds to the present results and the blue dotted line represents the literature results [22]. The right graph shows the same trajectory but compares it to rotations happening on the other wall.</figcaption></figure>

after the first interaction between the ellipsoid and the boundary of the computational
domain, at time $t \approx 0.45\text{s}$, are in good agreement.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the horizontal position of the center of the ellipsoid. The dotted line shows the literature results [23] and the red line represents the present results.</figcaption></figure>
## 5.2 Multiple objects

**Two disks falling in an incompressible fluid**
Collisions between two spherical bodies placed on a same vertical line and falling under the effect of gravity, in a channel filled with incompressible Newtonian fluid, are known as the *drafting, kissing and tumbling* phenomenon. We consider two disks of radius $r = 0.1\text{cm}$ and density $\rho_S = 1.01\frac{g}{\text{cm}^3}$ placed in a fluid domain of dimensions $[0\text{cm}, 2\text{cm}] \times [0\text{cm}, 8\text{cm}]$. To have a slight asymmetry, necessary to speed up the insurgence of the phenomenon, the coordinates of the upper disk $S1$ center are fixed at $(0.999\text{cm}, 7.2\text{cm})$ and those of the lower disk $S2$ center at $(1.0\text{cm}, 6.8\text{cm})$. The fluid density and viscosity are respectively given by $\rho_f = 1.0\frac{g}{\text{cm}^3}$ and $\mu = 0.01\frac{g}{\text{cm s}}$. The simulation is performed on the time interval $[0\text{s}, 5\text{s}]$ with a time step equal to $\Delta t = 0.001\text{s}$ and a mesh size set to $h = 0.01$. To model the *drafting, kissing and tumbling* phenomenon as in [24], one fixes the width of the collision zone at $\rho = 0.0225\text{cm}$ and uses stiffness parameters for the disk-disk and disk-fluid domain collisions given by $\epsilon = 7 * 10^{-5}$ and $\epsilon_F = 5 * 10^{-5}$. Figure 5 shows the results of the simulation. When both disks are next to each other, the motion of the lower disk reduces the fluid pressure behind it and therefore the resistance of the fluid for the upper disk. As a result, the upper body falls faster until it collides with the lower one. The collision forces cause the separation of the disks which then move in opposite directions. Figure 5 also compares the vertical and horizontal position of the two disks over time to literature [24]. It can be observed that the results of the collision phase have same behavior. In the present simulation, the bodies fall faster to the bottom, therefore the separation phase takes place sooner than in literature results. This explains the differences on the graphes.

**100 disks falling in an incompressible fluid**
We simulate the interactions between 100 circular particles in two dimensions. The particles of radius $r = 0.03125\text{cm}$ and density $\rho_S = 1.1\frac{g}{\text{cm}^3}$ are placed in a vertical channel of height 2cm and width 1cm which is filled with a fluid of density $\rho_f = 1.0\frac{g}{\text{cm}^3}$ and viscosity $\mu = 0.01\frac{g}{\text{cm s}}$. At the beginning of the simulation, the fluid and the particles are at rest. Then the particles move towards the bottom of the domain under the effect of gravity. The initial configuration is given by the first figure of 6. We run the simulation on a time interval $[0\text{s}, 5.5\text{s}]$ and use a time step fixed at $\Delta t = 0.002\text{s}$. The mesh size is given by $h = 0.01\text{cm}$. The collision force is defined on range $\rho = 0.025\text{cm}$. According to simulations of multiple objects interactions in the literature [25], the force intensity must be important to prevent the objects from overlapping. For this reason, we use stiffness parameters fixed at $\epsilon = 5 * 10^{-6}$ and $\epsilon_F = 10^{-7}$. Figure 6 shows the position of the 100 particles at four time instants. Due to the collision forces, the particles near the domain walls fall slower towards the bottom than the central particles. During the simulation the particles settle one onto the other on the bottom of the channel. At final time, almost all particles are in a stationary position and packed in a hexagonal lattice. Article [26] illustrates the same simulation but with different values of the parameters.

**Two objects in a flow**
For this application we consider the motion of two circular particles in a symmetric stenotic artery in two dimensions. In contrast with the previous test cases, the particles do not fall under the effect of gravity, but the motion of the
## 5.2 Multiple objects

**Two disks falling in an incompressible fluid**
Collisions between two spherical bodies placed on a same vertical line and falling under the effect of gravity, in a channel filled with incompressible Newtonian fluid, are known as the *drafting, kissing and tumbling* phenomenon. We consider two disks of radius $r = 0.1\text{cm}$ and density $\rho_S = 1.01\frac{g}{\text{cm}^3}$ placed in a fluid domain of dimensions $[0\text{cm}, 2\text{cm}] \times [0\text{cm}, 8\text{cm}]$. To have a slight asymmetry, necessary to speed up the insurgence of the phenomenon, the coordinates of the upper disk $S1$ center are fixed at $(0.999\text{cm}, 7.2\text{cm})$ and those of the lower disk $S2$ center at $(1.0\text{cm}, 6.8\text{cm})$. The fluid density and viscosity are respectively given by $\rho_f = 1.0\frac{g}{\text{cm}^3}$ and $\mu = 0.01\frac{g}{\text{cm s}}$. The simulation is performed on the time interval $[0\text{s}, 5\text{s}]$ with a time step equal to $\Delta t = 0.001\text{s}$ and a mesh size set to $h = 0.01$. To model the *drafting, kissing and tumbling* phenomenon as in [24], one fixes the width of the collision zone at $\rho = 0.0225\text{cm}$ and uses stiffness parameters for the disk-disk and disk-fluid domain collisions given by $\epsilon = 7 * 10^{-5}$ and $\epsilon_F = 5 * 10^{-5}$. Figure 5 shows the results of the simulation. When both disks are next to each other, the motion of the lower disk reduces the fluid pressure behind it and therefore the resistance of the fluid for the upper disk. As a result, the upper body falls faster until it collides with the lower one. The collision forces cause the separation of the disks which then move in opposite directions. Figure 5 also compares the vertical and horizontal position of the two disks over time to literature [24]. It can be observed that the results of the collision phase have same behavior. In the present simulation, the bodies fall faster to the bottom, therefore the separation phase takes place sooner than in literature results. This explains the differences on the graphes.

**100 disks falling in an incompressible fluid**
We simulate the interactions between 100 circular particles in two dimensions. The particles of radius $r = 0.03125\text{cm}$ and density $\rho_S = 1.1\frac{g}{\text{cm}^3}$ are placed in a vertical channel of height 2cm and width 1cm which is filled with a fluid of density $\rho_f = 1.0\frac{g}{\text{cm}^3}$ and viscosity $\mu = 0.01\frac{g}{\text{cm s}}$. At the beginning of the simulation, the fluid and the particles are at rest. Then the particles move towards the bottom of the domain under the effect of gravity. The initial configuration is given by the first figure of 6. We run the simulation on a time interval $[0\text{s}, 5.5\text{s}]$ and use a time step fixed at $\Delta t = 0.002\text{s}$. The mesh size is given by $h = 0.01\text{cm}$. The collision force is defined on range $\rho = 0.025\text{cm}$. According to simulations of multiple objects interactions in the literature [25], the force intensity must be important to prevent the objects from overlapping. For this reason, we use stiffness parameters fixed at $\epsilon = 5 * 10^{-6}$ and $\epsilon_F = 10^{-7}$. Figure 6 shows the position of the 100 particles at four time instants. Due to the collision forces, the particles near the domain walls fall slower towards the bottom than the central particles. During the simulation the particles settle one onto the other on the bottom of the channel. At final time, almost all particles are in a stationary position and packed in a hexagonal lattice. Article [26] illustrates the same simulation but with different values of the parameters.

**Two objects in a flow**
For this application we consider the motion of two circular particles in a symmetric stenotic artery in two dimensions. In contrast with the previous test cases, the particles do not fall under the effect of gravity, but the motion of the
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: The upper figures show the position of the two disks as well as the streamlines at four different time instants. The lower graphs compare the evolution of the disks horizontal and vertical position to literature [24]. The solid lines correspond to present results and the dotted one to literature results.</figcaption></figure>

particles is due to a pressure difference between the inlet and the outlet. This simulation is presented in [27] and [28].

The geometry, given in figure 7, is a channel of length $L = 32d$ and width $w = 8d$. The centerlines are represented by dotted lines. The diameter $d$ of the particles is set to $d = 8.5 \cdot 10^{-4} \text{cm}$. To create the stenosis, the authors of [27] and [28] add two symmetrical
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: The upper figures show the position of the two disks as well as the streamlines at four different time instants. The lower graphs compare the evolution of the disks horizontal and vertical position to literature [24]. The solid lines correspond to present results and the dotted one to literature results.</figcaption></figure>

particles is due to a pressure difference between the inlet and the outlet. This simulation is presented in [27] and [28].

The geometry, given in figure 7, is a channel of length $L = 32d$ and width $w = 8d$. The centerlines are represented by dotted lines. The diameter $d$ of the particles is set to $d = 8.5 \cdot 10^{-4} \text{cm}$. To create the stenosis, the authors of [27] and [28] add two symmetrical
<figure><img src="image_8.png" /><figcaption>Figure 6: Streamlines and position at four different time instants of 100 disks falling under the effect of gravity in an incompressible Newtonian fluid.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Geometry of symmetric stenotic artery</figcaption></figure>

protuberances shaped like semicircles. The radius of these protuberances depends on the width of the stenosis throat, which is fixed at $b = 1.75d$. The fluid viscosity is set to $\mu = 0.01 \text{ cm}^2 \text{s}$ and both fluid and particle densities are equal to $\rho_f = \rho_S = 1 \text{ g/cm}^3$. At time $t = 0$, the fluid and the particles are at rest. Then the particles move towards the protuberances due to the pressure difference $\Delta p = 541 \text{ Pa}$ between the inlet and outlet of the channel. The two particles are initially located at $8d$ to the left of the vertical centerline and asymmetrically with respect to the horizontal centerline: $h = 2d + \frac{d}{4000}$ and $h' = 2d$. This slight asymmetry allows them to cross the stenosis throat. The mesh size and time step of the simulation are respectively set to $h = 6.07 \times 10^{-5} \text{ cm}$
<figure><img src="image_8.png" /><figcaption>Figure 6: Streamlines and position at four different time instants of 100 disks falling under the effect of gravity in an incompressible Newtonian fluid.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Geometry of symmetric stenotic artery</figcaption></figure>

protuberances shaped like semicircles. The radius of these protuberances depends on the width of the stenosis throat, which is fixed at $b = 1.75d$. The fluid viscosity is set to $\mu = 0.01 \frac{\text{cm}^2}{\text{s}}$ and both fluid and particle densities are equal to $\rho_f = \rho_S = 1 \frac{\text{g}}{\text{cm}^3}$. At time $t = 0$, the fluid and the particles are at rest. Then the particles move towards the protuberances due to the pressure difference $\Delta p = 541\text{Pa}$ between the inlet and outlet of the channel. The two particles are initially located at $8d$ to the left of the vertical centerline and asymmetrically with respect to the horizontal centerline: $h = 2d + \frac{d}{4000}$ and $h' = 2d$. This slight asymmetry allows them to cross the stenosis throat. The mesh size and time step of the simulation are respectively set to $h = 6.07 * 10^{-5}\text{cm}$
and $\Delta t = 10^{-5}$s. The time interval is [0.0s, 0.01s]. Collision forces are applied when the particles are close to the protuberances. The width of the collision zone is fixed at $\rho = 0.0003 \text{cm}$ and the stiffness parameters of the particle-particle and particle-domain interactions are given by $\epsilon = 10^{-11}$ and $\epsilon_F = 7.5 * 10^{-12}$. Figure 8 shows the trajectory of the two particles: they have a similar motion (snapshots 1-9) until they get close to the stenosis throat (snapshot 10). Since the width of the stenosis throat is too small to allow both particles to pass simultaneously, the upper particle stops and changes direction due to collision forces. Snapshot 11 shows that it moves back to allow the lower particle passing the throat. Finally, the upper particle follows the lower particle through the artery (snapshots 11-17). The same particle behavior is present in the cited articles. The particle trajectory after passing the stenosis throat depends on the initial configuration and the definition of collision parameters.

<figure><img src="image_10.png" /><figcaption>Figure 8: Position of two particles at different time instants.</figcaption></figure>

### Three-sphere swimmer in a Stokes flow
For the last test case, we consider the tilted three-sphere swimmer placed close to the boundary of a horizontal channel in a Stokes regime. The radius of the swimmer spheres is set to $r = 1 \text{cm}$ and the length of its rods to $l = 10 \text{cm}$. The height and width of the computational domain are fixed at $h = 150 \text{cm}$ and $w = 40 \text{cm}$. The mass center of the central sphere is initially located at (15cm, 10cm), and the density of the swimmer is given by $\rho_S = 0.1 \frac{g}{\text{cm}^3}$. The fluid has a density equal to $\rho_f = 1 \frac{g}{\text{cm}^3}$ and a viscosity of $\mu = 1 \frac{\text{cm}^2}{\text{s}}$. The mesh size is set to
and $\Delta t = 10^{-5}$s. The time interval is [0.0s, 0.01s]. Collision forces are applied when the particles are close to the protuberances. The width of the collision zone is fixed at $\rho = 0.0003 \text{cm}$ and the stiffness parameters of the particle-particle and particle-domain interactions are given by $\epsilon = 10^{-11}$ and $\epsilon_F = 7.5 * 10^{-12}$. Figure 8 shows the trajectory of the two particles: they have a similar motion (snapshots 1-9) until they get close to the stenosis throat (snapshot 10). Since the width of the stenosis throat is too small to allow both particles to pass simultaneously, the upper particle stops and changes direction due to collision forces. Snapshot 11 shows that it moves back to allow the lower particle passing the throat. Finally, the upper particle follows the lower particle through the artery (snapshots 11-17). The same particle behavior is present in the cited articles. The particle trajectory after passing the stenosis throat depends on the initial configuration and the definition of collision parameters.

<figure><img src="image_10.png" /><figcaption>Figure 8: Position of two particles at different time instants.</figcaption></figure>

### Three-sphere swimmer in a Stokes flow
For the last test case, we consider the tilted three-sphere swimmer placed close to the boundary of a horizontal channel in a Stokes regime. The radius of the swimmer spheres is set to $r = 1 \text{cm}$ and the length of its rods to $l = 10 \text{cm}$. The height and width of the computational domain are fixed at $h = 150 \text{cm}$ and $w = 40 \text{cm}$. The mass center of the central sphere is initially located at (15cm, 10cm), and the density of the swimmer is given by $\rho_S = 0.1 \frac{g}{\text{cm}^3}$. The fluid has a density equal to $\rho_f = 1 \frac{g}{\text{cm}^3}$ and a viscosity of $\mu = 1 \frac{\text{cm}^2}{\text{s}}$. The mesh size is set to
$h = 0.3\text{cm}$ and the time step to $\Delta t = 0.1\text{s}$. Due to the orientation of the swimmer at $t = 0\text{s}$ and due to its swimming strategy, leading to a straight motion, the swimmer gets closer to the boundary of the channel. Once its right sphere arrives in the collision zone of width $\rho = 0.03\text{cm}$, lubrication forces, defined by a stiffness parameter set to $\epsilon_F = 5 \cdot 10^{-8}$, are applied onto the swimmer, who starts to change direction. The swimmer rotates until its left sphere enters the collision zone, where the repulsive force applied on this sphere forces the swimmer to move upwards, away from the boundary. This behavior of the three-sphere swimmer close to the boundary is shown in figure 9.

<figure><img src="image_11.png" /><figcaption>Figure 9: Behavior of the three-sphere swimmer close to the boundary of the computational domain.</figcaption></figure>
$h = 0.3\text{cm}$ and the time step to $\Delta t = 0.1\text{s}$. Due to the orientation of the swimmer at $t = 0\text{s}$ and due to its swimming strategy, leading to a straight motion, the swimmer gets closer to the boundary of the channel. Once its right sphere arrives in the collision zone of width $\rho = 0.03\text{cm}$, lubrication forces, defined by a stiffness parameter set to $\epsilon_F = 5 \cdot 10^{-8}$, are applied onto the swimmer, who starts to change direction. The swimmer rotates until its left sphere enters the collision zone, where the repulsive force applied on this sphere forces the swimmer to move upwards, away from the boundary. This behavior of the three-sphere swimmer close to the boundary is shown in figure 9.

<figure><img src="image_11.png" /><figcaption>Figure 9: Behavior of the three-sphere swimmer close to the boundary of the computational domain.</figcaption></figure>
# Reproducibility

The validation benchmarks and all applications illustrated in this paper are available in a public GitHub repository [29]. All the results can be reproduced.

# 6 Conclusions

In this paper, we presented a repulsive collision-avoidance model for arbitrary-shaped bodies, where the definitions of force direction and magnitude are based on the computation of the distance function using the narrow-band fast marching method. The model is validated against literature results, with which we find good agreement. Various applications of two and three-dimensional cases are illustrated. In the appendix, we present a tool for the insertion of solid bodies into meshed domains, which allows simulations where the conforming treatment of the fluid-solid interface is required. Determining the optimal values for the collision parameters is still an open problem which could be tackled by anticipating the moments of collision and controlling the displacement of the bodies. The extension of this work to elastic bodies and collisions with multiple contact points will be the subject of an upcoming paper.

# A Insertion of arbitrary solids into discretized fluid domains

In biological processes, such as particle transport in blood vessels, solid bodies move in geometrically complex domains. The reconstruction of such environments for numerical simulations is often image-based rather than relying on a CAD description, which complicates the mesh-fitted definition of solid bodies inside these environments. Hence, we have developed a tool allowing the insertion of arbitrary bodies in meshed domains with the identification of their geometric and material properties. The tool relies on the MMG library [30] feature that allows remeshing from a level-set function: given the initial triangulated environment and a level-set description of the solid body, it outputs a new triangulation that includes the solid, and where the interface between the two is conforming. In addition to this, we are able to insert multiple objects simultaneously, with prescribed position and orientation, and ensure that physical markers associated to the fluid and the solids are correctly transferred onto the new mesh.

The capabilities of the tool are illustrated in figures 10 and 11. Figure 10 shows a two-dimensional environment of complex geometry where a collection of objects have been successfully inserted without changes to their shapes. In particular all corners of the rectangular body have been recovered. Figure 11 illustrates a simpler three-dimensional domain where multiple bodies have been again successfully inserted. All three resulting meshes can be used for the simulation of fluid-solid interaction with collision treatment, as described in the first part of the paper.
# Reproducibility

The validation benchmarks and all applications illustrated in this paper are available in a public GitHub repository [29]. All the results can be reproduced.

# 6 Conclusions

In this paper, we presented a repulsive collision-avoidance model for arbitrary-shaped bodies, where the definitions of force direction and magnitude are based on the computation of the distance function using the narrow-band fast marching method. The model is validated against literature results, with which we find good agreement. Various applications of two and three-dimensional cases are illustrated. In the appendix, we present a tool for the insertion of solid bodies into meshed domains, which allows simulations where the conforming treatment of the fluid-solid interface is required. Determining the optimal values for the collision parameters is still an open problem which could be tackled by anticipating the moments of collision and controlling the displacement of the bodies. The extension of this work to elastic bodies and collisions with multiple contact points will be the subject of an upcoming paper.

# A Insertion of arbitrary solids into discretized fluid domains

In biological processes, such as particle transport in blood vessels, solid bodies move in geometrically complex domains. The reconstruction of such environments for numerical simulations is often image-based rather than relying on a CAD description, which complicates the mesh-fitted definition of solid bodies inside these environments. Hence, we have developed a tool allowing the insertion of arbitrary bodies in meshed domains with the identification of their geometric and material properties. The tool relies on the MMG library [30] feature that allows remeshing from a level-set function: given the initial triangulated environment and a level-set description of the solid body, it outputs a new triangulation that includes the solid, and where the interface between the two is conforming. In addition to this, we are able to insert multiple objects simultaneously, with prescribed position and orientation, and ensure that physical markers associated to the fluid and the solids are correctly transferred onto the new mesh.

The capabilities of the tool are illustrated in figures 10 and 11. Figure 10 shows a two-dimensional environment of complex geometry where a collection of objects have been successfully inserted without changes to their shapes. In particular all corners of the rectangular body have been recovered. Figure 11 illustrates a simpler three-dimensional domain where multiple bodies have been again successfully inserted. All three resulting meshes can be used for the simulation of fluid-solid interaction with collision treatment, as described in the first part of the paper.
<figure><img src="image_12.png" /><figcaption>Figure 10: Creation of two-dimensional bodies inside a complex discretized environment.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Two cases showing the insertion of bodies in three dimensions.</figcaption></figure>

## Acknowledgements

This work of the Interdisciplinary Thematic Institute IRMIA++, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program. The authors acknowledge the financial support of the French Agence Nationale de la Recherche (grant ANR-21-CE45-0013 project NEMO), and Cemosis. In addition, the authors thank A. Froehly and L. Cirrottola for their helpful comments and discussions.

## References

[1] R. Glowinski, Finite element methods for incompressible viscous flow, Vol. 9 of Handbook of Numerical Analysis, Elsevier, 2003.
<figure><img src="image_12.png" /><figcaption>Figure 10: Creation of two-dimensional bodies inside a complex discretized environment.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Two cases showing the insertion of bodies in three dimensions.</figcaption></figure>

## Acknowledgements

This work of the Interdisciplinary Thematic Institute IRMIA++, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program. The authors acknowledge the financial support of the French Agence Nationale de la Recherche (grant ANR-21-CE45-0013 project NEMO), and Cemosis. In addition, the authors thank A. Froehly and L. Cirrottola for their helpful comments and discussions.

## References

[1] R. Glowinski, Finite element methods for incompressible viscous flow, Vol. 9 of Handbook of Numerical Analysis, Elsevier, 2003.
[2] A. Wachs, M. Uhlmann, J. Derksen, D. P. Huet, Modeling of short-range interactions between both spherical and non-spherical rigid particles, in: Modeling Approaches and Computational Methods for Particle-Laden Turbulent Flows, Elsevier, 2023, pp. 217–264. doi:10.1016/B978-0-32-390133-8.00019-0.
URL https://linkinghub.elsevier.com/retrieve/pii/B9780323901338000190

[3] R. Jain, S. Tschisgale, J. Fröhlich, A collision model for DNS with ellipsoidal particles in viscous fluid, International Journal of Multiphase Flow 120 (2019) 103087. doi: 10.1016/j.ijmultiphaseflow.2019.103087.
URL http://arxiv.org/abs/2201.04456

[4] P. Singh, T. Hesla, D. Joseph, Distributed Lagrange multiplier method for partic- ulate flows with collisions, International Journal of Multiphase Flow 29 (3) (2003) 495–509. doi:10.1016/S0301-9322(02)00164-7.
URL https://linkinghub.elsevier.com/retrieve/pii/S0301932202001647

[5] A. Broms, A.-K. Tornberg, A Barrier Method for Contact Avoiding Particles in Stokes Flow, preprint (Jan. 2023).
URL http://arxiv.org/abs/2301.01666

[6] A. Lefebvre, Numerical simulation of gluey particles, ESAIM: Mathematical Modelling and Numerical Analysis 43 (1) (2009) 53–80. doi:10.1051/m2an/2008042.
URL http://www.esaim-m2an.org/10.1051/m2an/2008042

[7] M. N. Ardekani, P. Costa, W.-P. Breugem, L. Brandt, Numerical Study of the Sedimentation of Spheroidal Particles, International Journal of Multiphase Flow 87 (2016) 16–34. doi:10.1016/j.ijmultiphaseflow.2016.08.005.
URL http://arxiv.org/abs/1602.05769

[8] B. Muth, M.-K. Müller, M. Müller, P. Eberhard, S. Luding, Collision detection and administration methods for many particles with different sizes, in: 4th International Conference on Discrete Element Methods, DEM 4, Minerals Engineering Int., 2007, pp. 1–18.

[9] B. Vemuri, L. Chen, L. Vu-Quoc, X. Zhang, O. Walton, Efficient and Accurate Collision Detection for Granular Flow Simulation, Graphical Models and Image Processing 60 (6) (1998) 403–422. doi:10.1006/gmip.1998.0479.
URL https://linkinghub.elsevier.com/retrieve/pii/S1077316998904798

[10] M. Pouchol, A. Ahmad, B. Crespin, O. Terraz, A hierarchical hashing scheme for nearest neighbor search and broad-phase collision detection, Journal of Graphics, GPU, and Game Tools 14 (2) (2009) 45–59.

[11] C. Ericson, Real-time collision detection, Crc Press, 2004.

[12] J. A. Sethian, A fast marching level set method for monotonically advancing fronts., Proceedings of the National Academy of Sciences 93 (4) (1996) 1591–1595. doi: 10.1073/pnas.93.4.1591.
URL https://pnas.org/doi/full/10.1073/pnas.93.4.1591
[2] A. Wachs, M. Uhlmann, J. Derksen, D. P. Huet, Modeling of short-range interactions between both spherical and non-spherical rigid particles, in: Modeling Approaches and Computational Methods for Particle-Laden Turbulent Flows, Elsevier, 2023, pp. 217–264. doi:10.1016/B978-0-32-390133-8.00019-0.
URL https://linkinghub.elsevier.com/retrieve/pii/B9780323901338000190

[3] R. Jain, S. Tschisgale, J. Fröhlich, A collision model for DNS with ellipsoidal particles in viscous fluid, International Journal of Multiphase Flow 120 (2019) 103087. doi: 10.1016/j.ijmultiphaseflow.2019.103087.
URL http://arxiv.org/abs/2201.04456

[4] P. Singh, T. Hesla, D. Joseph, Distributed Lagrange multiplier method for partic- ulate flows with collisions, International Journal of Multiphase Flow 29 (3) (2003) 495–509. doi:10.1016/S0301-9322(02)00164-7.
URL https://linkinghub.elsevier.com/retrieve/pii/S0301932202001647

[5] A. Broms, A.-K. Tornberg, A Barrier Method for Contact Avoiding Particles in Stokes Flow, preprint (Jan. 2023).
URL http://arxiv.org/abs/2301.01666

[6] A. Lefebvre, Numerical simulation of gluey particles, ESAIM: Mathematical Modelling and Numerical Analysis 43 (1) (2009) 53–80. doi:10.1051/m2an/2008042.
URL http://www.esaim-m2an.org/10.1051/m2an/2008042

[7] M. N. Ardekani, P. Costa, W.-P. Breugem, L. Brandt, Numerical Study of the Sedimentation of Spheroidal Particles, International Journal of Multiphase Flow 87 (2016) 16–34. doi:10.1016/j.ijmultiphaseflow.2016.08.005.
URL http://arxiv.org/abs/1602.05769

[8] B. Muth, M.-K. Müller, M. Müller, P. Eberhard, S. Luding, Collision detection and administration methods for many particles with different sizes, in: 4th International Conference on Discrete Element Methods, DEM 4, Minerals Engineering Int., 2007, pp. 1–18.

[9] B. Vemuri, L. Chen, L. Vu-Quoc, X. Zhang, O. Walton, Efficient and Accurate Collision Detection for Granular Flow Simulation, Graphical Models and Image Processing 60 (6) (1998) 403–422. doi:10.1006/gmip.1998.0479.
URL https://linkinghub.elsevier.com/retrieve/pii/S1077316998904798

[10] M. Pouchol, A. Ahmad, B. Crespin, O. Terraz, A hierarchical hashing scheme for nearest neighbor search and broad-phase collision detection, Journal of Graphics, GPU, and Game Tools 14 (2) (2009) 45–59.

[11] C. Ericson, Real-time collision detection, Crc Press, 2004.

[12] J. A. Sethian, A fast marching level set method for monotonically advancing fronts., Proceedings of the National Academy of Sciences 93 (4) (1996) 1591–1595. doi: 10.1073/pnas.93.4.1591.
URL https://pnas.org/doi/full/10.1073/pnas.93.4.1591
[13] L. Berti, V. Chabannes, L. Giraldi, C. Prud'homme, Fluid-rigid body interaction using the finite element method and ALE formulation: framework, implementation and benchmarkingPreprint.

[14] L. Formaggia, A. Quarteroni (Eds.), Cardiovascular mathematics: modeling and simulation of the circulatory system, no. 1 in MS & A : modeling, simulation & applications, Springer, Milano, 2009.

[15] B. Maury, Direct Simulations of 2D Fluid-Particle Flows in Biperiodic Domains, Journal of Computational Physics 156 (2) (1999) 325-351. doi:10.1006/jcph.1999.6365.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999199963659

[16] H. Kanchi, A. Masud, A 3D adaptive mesh moving scheme, International Journal for Numerical Methods in Fluids 54 (6-8) (2007) 923-944. doi:10.1002/fld.1512.
URL https://onlinelibrary.wiley.com/doi/10.1002/fld.1512

[17] A. Najafi, R. Golestanian, A Simplest Swimmer at Low Reynolds Number: Three Linked Spheres, Physical Review E 69 (6) (2004) 062901. doi:10.1103/PhysRevE.69.062901.
URL http://arxiv.org/abs/cond-mat/0402070

[18] D. Wan, S. Turek, Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method, International Journal for Numerical Methods in Fluids 51 (5) (2006) 531-566. doi:10.1002/fld.1129.
URL https://onlinelibrary.wiley.com/doi/10.1002/fld.1129

[19] L. Wang, Z. Guo, J. Mi, Drafting, kissing and tumbling process of two particles with different sizes, Computers & Fluids 96 (2014) 20-34. doi:10.1016/j.compfluid.2014.03.005.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045793014001042

[20] K. Usman, K. Walayat, R. Mahmood, N. Kousar, Analysis of solid particles falling down and interacting in a channel with sedimentation using fictitious boundary method, AIP Advances 8 (6) (2018) 065201. doi:10.1063/1.5035163.
URL http://aip.scitation.org/doi/10.1063/1.5035163

[21] E. Zamora, L. Battaglia, M. Storti, M. Cruchaga, R. Ortega, Numerical and experimental study of the motion of a sphere in a communicating vessel system subject to sloshing, Physics of Fluids 31 (8) (2019) 087106. doi:10.1063/1.5098999.
URL http://aip.scitation.org/doi/10.1063/1.5098999

[22] Z. Xia, K. W. Connington, S. Rapaka, P. Yue, J. J. Feng, S. Chen, Flow patterns in the sedimentation of an elliptical particle, Journal of Fluid Mechanics 625 (2009) 249-272. doi:10.1017/S0022112008005521.
URL https://www.cambridge.org/core/product/identifier/S0022112008005521/type/journal_article
[13] L. Berti, V. Chabannes, L. Giraldi, C. Prud'homme, Fluid-rigid body interaction using the finite element method and ALE formulation: framework, implementation and benchmarkingPreprint.

[14] L. Formaggia, A. Quarteroni (Eds.), Cardiovascular mathematics: modeling and simulation of the circulatory system, no. 1 in MS & A : modeling, simulation & applications, Springer, Milano, 2009.

[15] B. Maury, Direct Simulations of 2D Fluid-Particle Flows in Biperiodic Domains, Journal of Computational Physics 156 (2) (1999) 325-351. doi:10.1006/jcph.1999.6365.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999199963659

[16] H. Kanchi, A. Masud, A 3D adaptive mesh moving scheme, International Journal for Numerical Methods in Fluids 54 (6-8) (2007) 923-944. doi:10.1002/fld.1512.
URL https://onlinelibrary.wiley.com/doi/10.1002/fld.1512

[17] A. Najafi, R. Golestanian, A Simplest Swimmer at Low Reynolds Number: Three Linked Spheres, Physical Review E 69 (6) (2004) 062901. doi:10.1103/PhysRevE.69.062901.
URL http://arxiv.org/abs/cond-mat/0402070

[18] D. Wan, S. Turek, Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method, International Journal for Numerical Methods in Fluids 51 (5) (2006) 531-566. doi:10.1002/fld.1129.
URL https://onlinelibrary.wiley.com/doi/10.1002/fld.1129

[19] L. Wang, Z. Guo, J. Mi, Drafting, kissing and tumbling process of two particles with different sizes, Computers & Fluids 96 (2014) 20-34. doi:10.1016/j.compfluid.2014.03.005.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045793014001042

[20] K. Usman, K. Walayat, R. Mahmood, N. Kousar, Analysis of solid particles falling down and interacting in a channel with sedimentation using fictitious boundary method, AIP Advances 8 (6) (2018) 065201. doi:10.1063/1.5035163.
URL http://aip.scitation.org/doi/10.1063/1.5035163

[21] E. Zamora, L. Battaglia, M. Storti, M. Cruchaga, R. Ortega, Numerical and experimental study of the motion of a sphere in a communicating vessel system subject to sloshing, Physics of Fluids 31 (8) (2019) 087106. doi:10.1063/1.5098999.
URL http://aip.scitation.org/doi/10.1063/1.5098999

[22] Z. Xia, K. W. Connington, S. Rapaka, P. Yue, J. J. Feng, S. Chen, Flow patterns in the sedimentation of an elliptical particle, Journal of Fluid Mechanics 625 (2009) 249-272. doi:10.1017/S0022112008005521.
URL https://www.cambridge.org/core/product/identifier/S0022112008005521/type/journal_article
[23] T.-W. Pan, R. Glowinski, G. P. Galdi, Direct simulation of the motion of a settling ellipsoid in Newtonian fluid, Journal of Computational and Applied Mathematics 149 (1) (2002) 71–82.

[24] Z.-G. Feng, E. E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems, Journal of Computational Physics 195 (2) (2004) 602–628. doi:10.1016/j.jcp.2003.10.013.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999103005758

[25] T.-W. Pan, D. D. Joseph, R. Bai, R. Glowinski, V. Sarin, Fluidization of 1204 spheres: simulation and experiment, Journal of Fluid Mechanics 451 (2002) 169–191.
doi:10.1017/S0022112001006474.
URL https://www.cambridge.org/core/product/identifier/S0022112001006474/type/journal_article

[26] L. H. Juárez, R. Glowinski, T. Pan, Numerical simulation of fluid flow with moving and free boundaries, SeMA Journal: Boletín de la Sociedad Española de Matemática Aplicada (30) (2004) 49–102.

[27] J. W. a. C. Shu, Particulate Flow Simulation Via a Boundary Condition-enforced Immersed Boundary-lattice Boltzmann Scheme, Communications in Computational Physics 7 (4) (2010) 793–812. doi:10.4208/cicp.2009.09.054.
URL http://global-sci.org/intro/article_detail/cicp/7655.html

[28] H. Li, H. Fang, Z. Lin, S. Xu, S. Chen, Lattice Boltzmann simulation on particle suspensions in a two-dimensional symmetric stenotic artery, Physical Review E 69 (3) (2004) 031919. doi:10.1103/PhysRevE.69.031919.
URL https://link.aps.org/doi/10.1103/PhysRevE.69.031919

[29] C. Prud'homme, V. Chabannes, T. Metivet, R. Hild, Trophime, S. Abdoulaye, T. Saigre, P. Ricka, L. Berti, C. Van Landeghem, feelpp/feelpp. doi:10.5281/ZENODO.5718297.

[30] G. Balarac, F. Basile, P. Bénard, F. Bordeu, J.-B. Chapelier, L. Cirrottola, G. Caumon, C. Dapogny, P. Frey, A. Froehly, G. Ghigliotti, R. Laraufie, G. Lartigue, C. Legentil, R. Mercier, V. Moureau, C. Nardoni, S. Pertant, M. Zakari, Tetrahedral remeshing in the context of large-scale numerical simulation and high performance computing, MathematicS In Action 11 (1) (2022) 129–164. doi:10.5802/msia.22.
URL https://msia.centre-mersenne.org/articles/10.5802/msia.22/
[23] T.-W. Pan, R. Glowinski, G. P. Galdi, Direct simulation of the motion of a settling ellipsoid in Newtonian fluid, Journal of Computational and Applied Mathematics 149 (1) (2002) 71-82.

[24] Z.-G. Feng, E. E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems, Journal of Computational Physics 195 (2) (2004) 602-628. doi:10.1016/j.jcp.2003.10.013.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999103005758

[25] T.-W. Pan, D. D. Joseph, R. Bai, R. Glowinski, V. Sarin, Fluidization of 1204 spheres: simulation and experiment, Journal of Fluid Mechanics 451 (2002) 169-191.
doi:10.1017/S0022112001006474.
URL https://www.cambridge.org/core/product/identifier/S0022112001006474/type/journal_article

[26] L. H. Juárez, R. Glowinski, T. Pan, Numerical simulation of fluid flow with moving and free boundaries, SeMA Journal: Boletín de la Sociedad Española de Matemática Aplicada (30) (2004) 49-102.

[27] J. W. a. C. Shu, Particulate Flow Simulation Via a Boundary Condition-enforced Immersed Boundary-lattice Boltzmann Scheme, Communications in Computational Physics 7 (4) (2010) 793-812. doi:10.4208/cicp.2009.09.054.
URL http://global-sci.org/intro/article_detail/cicp/7655.html

[28] H. Li, H. Fang, Z. Lin, S. Xu, S. Chen, Lattice Boltzmann simulation on particle suspensions in a two-dimensional symmetric stenotic artery, Physical Review E 69 (3) (2004) 031919. doi:10.1103/PhysRevE.69.031919.
URL https://link.aps.org/doi/10.1103/PhysRevE.69.031919

[29] C. Prud'homme, V. Chabannes, T. Metivet, R. Hild, Trophime, S. Abdoulaye, T. Saigre, P. Ricka, L. Berti, C. Van Landeghem, feelpp/feelpp. doi:10.5281/ZENODO.5718297.

[30] G. Balarac, F. Basile, P. Bénard, F. Bordeu, J.-B. Chapelier, L. Cirrottola, G. Caumon, C. Dapogny, P. Frey, A. Froehly, G. Ghigliotti, R. Laraufie, G. Lartigue, C. Legentil, R. Mercier, V. Moureau, C. Nardoni, S. Pertant, M. Zakari, Tetrahedral remeshing in the context of large-scale numerical simulation and high performance computing, MathematicS In Action 11 (1) (2022) 129-164. doi:10.5802/msia.22.
URL https://msia.centre-mersenne.org/articles/10.5802/msia.22/