An immersed method based on cut-cells for the simulation of 2D incompressible fluid flows past solid structures

TL;DR

This paper introduces a cut-cell method for simulating 2D incompressible flows around solid obstacles, combining finite volume discretization with a MAC scheme on Cartesian grids, demonstrating efficiency at various Reynolds numbers.

Contribution

The paper presents a novel immersed boundary method using cut-cells and MAC scheme for efficient 2D incompressible flow simulation around static and moving obstacles.

Findings

Effective simulation of flow around a circular cylinder at Re=1000 and 3000.

Successful application to flow around a moving rigid body at Re=800.

Validation of the method's efficiency and accuracy through numerical results.

Abstract

We present a cut-cell method for the simulation of 2D incompressible flows past obstacles. It consists in using the MAC scheme on cartesian grids and imposing Dirchlet boundary conditions for the velocity field on the boundary of solid structures following the Shortley-Weller formulation. In order to ensure local conservation properties, viscous and convecting terms are discretized in a finite volume way. The scheme is second order implicit in time for the linear part, the linear systems are solved by the use of the capacitance matrix method for non-moving obstacles. Numerical results of flows around an impulsively started circular cylinder are presented which confirm the efficiency of the method, for Reynolds numbers 1000 and 3000. An example of flows around a moving rigid body at Reynolds number 800 is also shown, a solver using the PETSc-Library has been prefered in this context to solve the linear systems.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}-\frac{1}{Re}\Delta \mathbf{u}+\nabla \cdot (\mathbf{u}\otimes \mathbf{u})+\nabla p=0,$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}(\mathbf{x},t=0)={\mathbf{u}}_0,$$

$$\mathbf{u}=0\text{mbox}on\partial {\Omega }^F$$

$$\frac{3{\tilde{\mathbf{u}}}^{k+1}-4{\mathbf{u}}^k+{\mathbf{u}}^{k-1}}{2\delta t}-\frac{1}{Re}\Delta {\tilde{\mathbf{u}}}^{k+1}+\nabla P^k=-2\nabla \cdot ({\mathbf{u}}^k\otimes {\mathbf{u}}^k)+\nabla \cdot ({\mathbf{u}}^{k-1}\otimes {\mathbf{u}}^{k-1})$$

$$\begin{array}{rl}& \frac{{\mathbf{u}}^{k+1}-{\tilde{\mathbf{u}}}^{k+1}}{\delta t}+\frac{2}{3}\nabla (P^{k+1}-P^k)=0,\\ & \nabla \cdot {\mathbf{u}}^{k+1}=0,({\mathbf{u}}^{k+1}-{\tilde{\mathbf{u}}}^{k+1})\cdot \mathbf{n}=0\text{on }\Gamma .\end{array}$$

$$AX=Z$$

$$GX=Z-QY$$

$$QY=(A-G)X.$$

$$\left(I_{n_c}+M{G}^{-1}Q\right)Y=M{G}^{-1}Z$$

$$\frac{\partial \mathbf{u}}{\partial t}+(1,0)\cdot \nabla \mathbf{u}=0$$

$$\mathbf{F}_{b}\,=\,\int_{\Gamma_{S}}\Bigl{(}-p\mathbf{n}+\frac{1}{\textrm{Re}}\frac{\partial\mathbf{u}}{\partial n}\Bigr{)}ds.$$

$$C_{\text{pred}}=4\sqrt{\frac{2\pi }{\text{Re}T}}+\frac{2\pi }{\text{Re}}\left(9-\frac{15}{\sqrt{\pi }}\right).$$

$${\mathbf{u}}_{\text{body}}(t)=\left(2\sin \left(\frac{t}{2}\right);0\right)$$

Full text

1. Introduction

For some decades, many researchers and engineers have been considering the numerical solution of fluid flows, for different kind of fluids and different geometries. With the increasing performance of super-computers, it has been possible to tackle more and more challenging problems, for higher Reynolds numbers and complex geometries. Several different discretization techniques can be used to consider these problems: Finite Element Methods, Finite Volumes Methods, Spectral Methods, and Finite Difference Methods. The MAC scheme on cartesian grids [15] can be viewed both as a Finite Volume Method or a Finite Difference Method on staggered grids, and is adapted to 2D or 3D flows in simple geometries, for example for lid-driven cavity or backward facing step. To take into account some obstacles in the flows (or complex geometries), immersed boundary techniques have been developped by Peskin in the 80’s ([26, 27]), consisting in using Dirac functions to model the interacting force between the fluid and the solid structure. These methods have inspired many authors in the following years, Mohd-Yusof has combined them with the use of B-Splines ([24]) in his momentum forcing methods to consider complex geometries. The main advantage of these techniques is that the forcing term does not change the spatial operators, making them quite easy to implement (see [23] for a review, and refences therein). As an alternative, Bruno et al. have developped penalization techniques to inforce suitable boundary conditions [1]. Similar techniques have also been investigated by Maury et al. [17] and justified from a mathematical point of view in [21]. These methods have been shown to be efficient in the context of several particles in a flow [19], and when considering possible collisions between them [32].

Arbitrary Lagrangian Eulerian (ALE) methods have been developped for flows in geometries which vary in time (see [28, 29, 33] where authors use some of the ideas of [5]). The aim is to formulate the equation in a fixed reference domain, by using a mapping from the reference domain $\Omega(0)$ to the domain $\Omega(t)$ occupied by the fluid at time $t$. The position of the moving bodies, which correspond to the boundary of $\Omega(t)$, being available, the velocity field of these bodies defined on $\partial\Omega(t)$ have to be extended to $\Omega(t)$. Once this is done (generally with harmonic extensions), the equations are written in the reference domain by using the chain-rule formula.

For problems involving non-rigid bodies, Roshchenko et al. (see [30]) have used splitting methods to solve first the evolution of the velocity field in the fluid, and then to consider the deformation of the body. These ideas of splitting the model can be viewed as similar to the projection techniques (see [11], and [14] for a review and references therein).

The method presented in this work joints another family of methods, called cut-cell methods. The idea of these methods is to modify the discretization of the Navier-Stokes Equations in the cells cut by the immersed boundary (see [34, 22, 31, 12]). One can discretize the equations on smallest cells obtained by intersecting the grid-cell with the domain occupied by the fluid, or one can merge these smallest cells with a neighbouring one. These methods can be combined with the levelset methods to track the boundary of the fixed or moving body (see [25]). These ideas are used in the present work: the body in the fluid is represented by a levelset function, and the location of the velocity components are modified in the cut-cells (see [7]), the pressure remaining placed at the center of cartesian grid cells for both fluid-cells and cut-cells. For the Laplacian of the velocity, the classical five-point approximation must be replaced by a local 6-point formula, for which the truncation error is only first order. But as in [20], global second order convergence of the method is recovered. This second-order convergence for the velocity and the pressure with our cut-cell scheme has been obtained for the flow past a circular cylinder at Re=40 in [16] by comparing with the reference solution proposed in [13].

The paper is organized as follows: Section 2 is devoted to the presentation of the problem. The Navier-Stokes Equations are considered in a 2D geometry, which is supposed to be fixed for the sake of clarity. We also introduce there the notation for the grids, and detail space discretization. In Section 3, we give some information about computational aspects. In the case of fixed domain, a fast solver adapted from the capacitance matrix method (see [10, 9]) is used to solve the linear systems for both components and for the pressure. We also show that the method can be adapted in the case of moving domain. In this context, the preprocessing step of the capacitance matrix method would need to be done at each time-iteration which would then increase the CPU time. Therefore, we have prefered for this case a parallel version of an algebraic multigrid method (HYPRE BoomerAMG) implemented using the PETSc Fortran library (see [2, 3]).

Section 4 is then devoted to numerical tests, where we compare our results with theoretical predictions and other numerical results available in the literature. A conclusion is given in section 5.

2. The Setting of the Problem

2.1. Flows past obstacles

We consider a flow in a two-dimensional domain $\Omega=(0,L)\times(0,H)$ which contains a domain $\Omega^{S}$ occupied by the solid which is supposed to be fixed for the sake of simplicity. We denote then $\Omega^{F}=\Omega\setminus\Omega^{S}$ the domain occupied by the fluid (see Fig. 1).

The velocity field in the fluid satisfies the Navier-Stokes equations, with no-slip boundary conditions. We consider then the problem:

[TABLE]

where ${\mathbf{u}}({\mathbf{x}},t)=(u,v)$ is the velocity field at ${\mathbf{x}}=(x,y)\in\Omega^{F}$ at time $t>0$, ${\mathbf{u}}_{0}$ is the initial condition and $Re$ is the Reynolds number. We impose homogeneous Dirichlet boundary conditions for the velocity field on $\partial{\Omega^{F}}$:

[TABLE]

We mention that non-homogenous Dirichlet boundary conditions can also be treated with the method presented here.

2.2. Discretization

For the time-discretization of (1)–(3), we use a second-order backward difference (BDF2) projection scheme. In a first step, the velocity field is advanced in time with a semi-implicit scheme decoupling the velocity and pressure unknowns. Then, the intermediate velocity is projected in order to obtain a free-divergence velocity field.

Let $\delta t>0$ stand for the time step and $t^{k}=k\,\delta t$ discrete time values. Let us consider that $(\mathbf{u}^{j},\,P^{j})$ are known for $j\leq k$. The computation of $(\mathbf{u}^{k+1},\,P^{k+1})$ needs two steps:

[TABLE]

with homogeneous Dirichlet boundary condition for $\widetilde{\mathbf{u}}^{k+1}$.

Then the intermediate velocity field $\widetilde{\mathbf{u}}^{k+1}$ is projected in the free-divergence space to get $\mathbf{u}^{k+1}$:

[TABLE]

For the spatial discretization, we modify the MAC scheme near the boundary by changing the location of the unknowns of the velocity components for the cells cut by the solid as depicted on Fig. 2 , the pressure unknowns remaining in their original place (see [7] for more details).

To discretize the Laplacian in (5), we must replace the five-point formula by a six-point discretization. For the convective terms, the fluxes are computed at the midle of the vertical and horizontal edges (see Fig. 3).

For the pressure, linear interpolation are used rather than changing the location of the pressure unknowns. The same kind of linear interpolation is used to get consistant evaluation of the pressure gradient in (6).

Although the truncation error is only first order in space for the resulting numerical scheme, the second order accuracy is recovered which is due to a superconvergence phenomena analoguous to those proven in [20]. This second order has been observed in [16] by showing results in comparison with those of [13].

3. Computational Aspects

3.1. A fast parallel direct solver to treat fixed solid strutures

When considering fixed solid bodies, the use of a direct solver with a preprocessing procedure is efficient. Once the preprocessing computations have been done, the cost paid to solve the linear systems is about twice the case of a numerical simulation in the same computational domain without obstacles. We summarize hereafter the fast direct solver derived from the capacitance matrix method and adapted to the case of non uniform grids (see [7] and [8] for the details) which has been implemented in our code.

After spatial discretization of the Navier-Stokes equations, one linear equation is obtained per node in the part of the computational domain filled by the fluid and per unknown that is $u,\,v$ and $p$. We complete these sets of linear equations by adding similar ones for nodes of the cartesian grid lying inside the solid obstacle but with zero as right-hand side. The unknowns corresponding to mesh points in $\Omega_{h}^{S}$ are fictitious ones. As in $\Omega^{F}$, the numerical scheme accounts for the boundary conditions on $\Gamma_{h}^{S}$, the fluid unknowns are independant to the solid ones. We therefore obtain linear algebraic systems defined on the whole cartesian grid with $n_{x}\times n_{y}$ mesh points whose sizes are $(n_{x}-1)\times n_{y}$ for $u$, $n_{x}\times(n_{y}-1)$ for $v$ and $(n_{x}-1)\times(n_{y}-1)$ for $p$. All three linear systems are similar in nature : the resulting matrices have similar structures with five or six non-zero coefficients per row.

Let us consider one of these linear system. We denote by $N$ its size and by $A\in\mathcal{M}_{N}(\mathbb{R})$ its matrix. Then at each time iteration, we have to solve a linear system

[TABLE]

with the right-hand side $Z$ computed from the velocity and the pressure at previous time steps. As it is mentioned above, the matrix $A$ is non-symmetric. Let us consider now the matrix $G$ obtained with the same discretization on the whole computational domain $\Omega$ totally filled by a fluid that is with no obstacles. The matrices $A$ and $G$ differ only on rows corresponding to computational meshes for which the five-point stencil interacts with a cut-cell. Let us denote by $n_{c}$ this total number of rows, namely rows such that $A-G$ have non-vanishing coefficients. The efficiency of our direct solver is due to the fact that $n_{c}$ is small compared with $N$ and that the non-zero coefficients on each row of $A-G$ is bounded. The linear system (7) can be rewritten as

[TABLE]

where $Q$ is a matrix of dimensions $N\times n_{c}$ with one non-vanishing coefficient per column, equal to one, and $Y\in\mathbb{R}^{n_{c}}$ such that

[TABLE]

It can be easily shown, using $Q^{t}Q=I_{n_{c}}$, that $Y$ is solution of the following linear system

[TABLE]

with $M=Q^{t}(A-G)$. The matrix $I_{n_{c}}\,+\,M\,G^{-1}Q$ is a non-singular matrix (see [8] for a proof) of size $n_{c}$.

Based on these relations, the algorithm implemented to solve (7) consists in a preprocessing step where the matrix $I_{n_{c}}\,+\,M\,G^{-1}Q$ is factorized (we use a $LU$-factorization) followed by

1. i)

Compute $Z$ and solve $GW=Z$ ; 2. ii)

Compute $MW$ and solve (9) ; 3. iii)

Compute $QY$ and solve $GX\,=\,Z\,-\,QY$.

Recalling that $G$ is the matrix corresponding to the standard MAC scheme on the whole computational mesh, steps i) and iii) can be performed by using any efficient solvers available on cartesian grids. In the present work, we use Discrete Fourier transforms in the vertical direction (where the mesh is uniform) combined with $LU$-factorizations of the resulting tridiagonal systems.

The parallel version of this direct solver is based on explicit communications performed by calling functions of the MPI library. The main feature of MPI is that a parallel application consists in running $p$ independant processes which may be executed on different computers, processors or cores. These processes can exchange datas by sending/receiving messages via a network connecting all the involved computing units.

The first step when developping a parallel algorithm is to define a suitable and efficient splitting of the datas among the MPI processes : each MPI process will treat datas associated with a part of the total computational mesh. For our problem, this choice is straightforward and is related to the algorithm used to solve the linear systems. Indeed, it is much easier to implement a parallel resolution of tridiagonal linear systems rather than a parallel version of the DFT. Therefore, the parallel version of the code is based on a splitting of the datas along the horizontal axis, so that each MPI process works with a vertical slice of the computational mesh as it is illustrated on Figure 4.

In the framework of finite volume or finite difference schemes on cartesian grids, the explicit computation of spatial derivatives is local and involves very few communications. The only tricky part concerns the resolution of the linear systems. The step ii) of the direct solver described in the previous section consists in solving a linear system involving the matrix $I_{n_{c}}\,+\,M\,G^{-1}Q$. As the $LU$-factorization of this matrix has been computed and stored in a pre-processing step at the beginning of the time iterations, we have to solve two triangular systems which can not be efficiently performed on parallel computers. As $n_{c}$ is small compared to the size of the global problem, we choose to dedicate this task to one given MPI process (fixed in advance) per unknown, that is $u,v$ and $p.$ Once these linear systems are solved, the resulting vectors are scattered from these MPI processes to the other ones.

As it was previously mentioned, linear systems of steps i) and ii) are solved by first applying a DFT in the $y$-direction : these computations are independant and can be performed without any communications due to the distribution of datas among the MPI processes. This results in a collection, one per grid point in the $y$-direction, of independant tridiagonal linear systems connecting all nodes of the mesh in the $x$-direction. A parallel direct solver based on the divide and conquer approach (DAC) for tridiagonal matrices has been implemented (see [6]). The DAC method, applied to solve one tridiagonal linear system on $n_{p}>1$ MPI processes, consists in splitting the tridiagonal matrix into $n_{p}$ independant blocks (one per MPI process). The solutions of these systems have to be corrected in order to recover the solution of the global system. These corrections correspond to $2n_{p}-1$ values which are solutions of a tridiagonal linear system of size $2n_{p}-1$. This phase of the DAC method is sequential and has to be performed on one process inducing a useless waiting time for the other processes. However, as we have to solve $n_{y}$ such systems simultaneously, this sequential part can be distributed among all the $n_{p}$ processes. In this context, the DAC algorithm leads to an efficient parallel code.

The parallel code has a good level of performance : less than $15\%$ of the CPU time is spent in communications between MPI processes. The sequential part performed on one process represents a negligible amount of CPU time. The computations presented here have been performed on a DELL cluster using up to 32 cores of Xeon processors. A low latency bandwith network connects the cluster nodes.

3.2. Iterative solver for the case of moving bodies

For solid bodies moving in a computational domain filled by a fluid, as the case considered in Section 4.2, cut-cells may change at each time iteration. Therefore, the coefficients of the matrices of the linear systems for the velocity components, issued from the discretizetion of the momentum equations, and for the pressure increment computed in the projection step of the time scheme, have to be recomputed at each time step. In that context, the use of the fast direct solver described in the previous section is cumbersome and inefficient except on coarse meshes. In order to be able to treat such configurations, we have implemented a PETSc version [2, 3] of our cut-cell scheme. The main advantage of the PETSc Library is that, in a parallel programming environment based on MPI, many iterative solvers combined with different preconditionners can be used. The choice can be made at run time.

4. Numerical Results

4.1. Flow past a circular cylinder at $\textrm{Re}=1000$ and $3000$

In this section, we present numerical simulations, performed with the parallel direct solver method described in Section 3.1. We consider the case of flows past a fixed circular cylinder of diameter $D$. The Reynolds number is defined based on the diameter $D$ of the cylinder, i.e. $\textrm{Re}=U_{\infty}D/\nu$ where $U_{\infty}$ is the horizontal free stream velocity. As non-dimensional time, we consider $T=2U_{\infty}t/D$.

A circular cylinder of diameter equals to unity is centered at the origin of the computational domain $\Omega=(-L_{x},L_{x})\times(-L_{y},L_{y})$. As boundary conditions, a uniform velocity profile $\mathbf{u}(\mathbf{x}=-L_{x},t)=(1,0)$ is imposed at the inflow and a convective boundary condition is applied at the exit, namely the convective equation

[TABLE]

is solved at $x=L_{x}$. On the top and bottom boundaries, that is $y=\pm L_{y}$, slip boundary conditions are used that is $\frac{\partial u}{\partial n}=0$ and $v=0$. Finally, no-slip ($\mathbf{u}|_{\Gamma_{S}}=\mathbf{0}$) boundary condition is applied on the surface of the obstacle.

For this problem important quantities reflecting the dynamics of the vortices formed in the vicinity of the solid boundary and developping at the rear of the cylinder are the pressure drag and lift coefficients. They are derived from the total drag force on the body, which is computed as

[TABLE]

The pressure drag and lift coefficients $C_{p}$ and $C_{\ell}$ are given by $C_{p}=2\mathbf{F}_{b}\cdot\mathbf{e}_{x}$ and $C_{\ell}=2\mathbf{F}_{b}\cdot\mathbf{e}_{y}$ and the (total) drag coefficient is $C_{d}=C_{p}+C_{\ell}$. Starting with a flow at rest, the drag coefficient behaves as $T^{-1/2}$ in the early stage of the development of the vortices. This square-root singularity has been theoretically predicted by Bar-Lev and Yang in [4]. They have derived the following expression for the (total) drag coefficient

[TABLE]

As whown in Figure 5, the values obtained with the cut-cell scheme on a grid with $4096\times 8192$ mesh points discretizing the domain $\Omega=(-10,10)^{2}$ perfectly match the theoretical curve drawing (12) on the time interval $T\in[0,0.2]$ for $\textrm{Re}=1000$. Our cut-cell method captures the square-root singularity of the drag coefficient. This simulation has been run using $16$ MPI processes. On longer time interval, namely $T\in[0,5]$, the results are in good agreement with those obtained by Koumoutsakos and Leonard in [18] with a vortex method. In order to test the grid convergence of these results, the same simulation has been conducted on a grid with two times more points in both spatial directions, that is $8192\times 16384$ mesh points, in the same computational domain. In that case, $32$ MPI processes have been used. Both results are almost indistinguishable on Figure 6 indicating that the coarser resolution is enough to capture the essential features of the flow at this Reynolds number.

The development of the flow around the impulsively started cylinder at $\textrm{Re}=1000$ can be seen on Figure 7. In the early stage $T\leq 1$, a primary vortex develops in the vinicity of the boundary at the rear of the cylinder. Then for $T\in[1,2]$ a secondary vortex appears trying to move insight the primary vortex and to push it away from the solid boundary ($T\geq 4$). A tertiary vortex is visible at $T=3$ which remains sticked to the boundary constrained by the two other vortices having more strength. These results compare well with the same flow representations shown in [18]. As expected on short time interval ($T\leq 5$) the flow remains symmetric.

At $\textrm{Re}=3000$, the time evolution of the drag coefficient plotted on Figure 8 exhibits also the square-root singularity on short time interval and is in good agreement with the results of Koumoutsakos and Leonard. As expected, the drag coefficients remains almost constant for $T\in[2,3]$ (see [18]). Note that on this time interval, a small difference exists between the results computed on the two different grids. However, the coarser simulation is fine enough to capture the flow dynamics at $\textrm{Re}=3000$. The mesh size of the coarser grid, which is constant in the vicinity of the solid boundary, is $h=20/8192\approx 2.44\times 10^{-3}$. Note that the size of the computational domain and the boundary conditions imposed at the exit may influence the results. Estimating the values for $L_{x}$ and $L_{y}$ required so that the numerical results being of the order of the numerical scheme error, namely $\textrm{O}(h^{2})$, is an open question. This will be addressed in further works. At this Reynolds number, a scenario similar than that at $\textrm{Re}=1000$ can be observed on Figure 9 with the development of three vortices in the early stage of the flow dynamics. The secondary vortex penetrates further inside the primary vortex aera and the tertiary vortex has more strength as it could be expected with less effects of the viscous forces at $\textrm{Re}=3000$. Again, an overall good agreement is found with the vortcity contours shown in [18] at the same Reynolds number and time $T$.

As previously mentioned, the flow remains symmetric at the beginning of the simulations for these flows around an impulsively started cylinder. By carrying the time integration over a much longer time interval $T\in[0,200]$ instabilities due to round-off errors and to the nonlinearity of the system develop so that the flow becomes non symmetric for $T\geq 100$ at $\textrm{Re}=1000$ and $T\geq 50$ at $\textrm{Re}=3000$ as it can be seen on Figures 10 and 11 representing the time history of the drag coefficient. After a transient period, an increase of the drag coeffient is observed which stabilizes and oscillates around a mean value.

4.2. Flow around moving bodies

The purpose of this section is to show that the present numerical method is also able to simulate incompressible flows around moving bodies. Let us consider a cylinder which starts to move impulsively at $t=0$ with the sinusoidal translational motion

[TABLE]

in a fluid initially at rest for Reynolds number $800$. We suppose that the fluid is confined within a rectangular computational domain $\Omega=[-3;3]\times[-1;1]$ with no-slip boundary condition on $\partial\Omega^{F}$. The diameter of the cylinder is equal to $1$ and it is initially centered at the origin. The boundary condition (13) at the body surface $\partial\Omega^{S}$ is enforced through the non exhaustive following right hand side terms which vanish in the case of fixed obstacle : $u(\mathbf{\kappa}_{i,j}^{u,S})$, $u(\xi_{i,j-1/2}^{S},y_{j-1/2})$, $\mathbf{u}_{\textrm{body}}(\mathbf{\kappa}_{i,j}^{S})$ and so on. This terms are respectively taking into account in the convective terms, Laplacian operator and continuity equation. More details can be found in [7].

As the obstacle moves from one time step to another, we have to update matrices of the linear systems corresponding to Poisson and momentum equations at each iteration. Therefore, in such configuration, the direct solver requires much more CPU time compared to some iterative solver, which does not require a preprocessing step. For large problems, the faster solver we have found is an algebraic multigrid method (HYPRE BoomerAMG) implemented using the PETSc library [2, 3].

A constant mesh size $h=5\times 10^{-3}$ is used in both directions and the value of the time step, satisfying a CFL stability condition, is $10^{-3}$. As shown in Figure 12, vortices interact with each other and also with the boundaries. The flow remains perfectly symmetric until $t=6\pi$, thereafter the symmetry of the flow is lost due to rounding errors inherent in computer calculations (see Figure 13).

5. Conclusion

We have presented a cut-cell method for the numerical solution of flows past obstacles. We have detailed the numerical method and the computational aspects for fixed obstacles, and shown numerical results for fixed and moving rigid bodies. The parallel version of the algorithm presented here allows computations of flows at Reynolds number up to 3000. The numerical tests confirm some results of the literature, a good agreement is observed with the numerical simulations in [18]. We have also shown that the computation of the drag coefficient matches the theoretical square-root singularity predicted by [4]. The choice of the size of the box (compared with the grid size $h$) is one of the questions that we would like to investigate with this method in further works, we also would like to deal with rigid bodies following the fluid flow.

Acknowledgements

This research was partially supported by the French Government Laboratory of Excellence initiative n${}^{\textrm{o}}$ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Fund. The numerical simulations have been performed on a DELL cluster with $32$ processors Xeon E2650v2 ($8$ cores), $1$ To of total memory and an infiniband (FDR 56Gb/s) connecting network.