Michel Belliard
email: michel.belliard@cea.fr

M Chandesris

J Dumas

Y Gorsse

Dominique Jamet

C Josserand

B Mathieu

An analysis

Keywords: Immersed Boundary Method, Direct Forcing method, Moving body, Spurious force oscillations, Second order accurate numerical method, Fluid-structure interaction problems

Introduction

Fluid-structure interaction (FSI) problems are widely present in industrial applications as aircraft wings, heatexchanger tube bundles, bridges, off-shore platforms, . . . But they are tricky to simulate due to the complexity of the physical phenomena at stake and the presence of time varying geometries. For instance, the heat-exchanger tube bundles of the nuclear industry are subject to fluid-elastic coupling forces, which depend on the bundle geometry, the fluid nature (single-phase or two-phase) and on the fluid velocity. To manage numerically these complex problems, it would be tempting to use the well-known body-fitted approach. In this case, the boundary conditions (BC), that are critical for FSI issues, are exactly imposed on the fluid-structure interface. The Arbitrary Lagrangian Eulerian framework is a prototype of this class of methods, see for instance [START_REF] Farhat | Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution[END_REF][START_REF] Farhat | The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids[END_REF]. However accurate this approach may be, it is sometimes not possible to use it to handle complex industrial problems involving large motions and/or deformations of the bodies, eventually with topological changes, which require complex numerical schemes to deal with the re-meshing issue and particular geometries.

Another approach consists in using Fictitious Domain Methods [START_REF] Saul'ev | On the solution of some boundary value problems on high performance computers by fictitious domain method[END_REF]. The general idea is to consider only a fluid domain in which the solid boundaries are immersed. This can be done by the local modification of the computational scheme or by adding a supplementary term in the governing equations of the fluid. For the former class of methods, we can cite the Immersed Interface Method [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], the "Cut-Cell" or Cartesian Grid methods [START_REF] Yang | A cartesian cut cell method for compressible flows Part A: Static body problems[END_REF][START_REF] Johansen | A Cartesian grid embedded boundary method for Poisson's equation on irregular domains[END_REF] and the Jump Embedded Boundary Condition approach [START_REF] Ph | A model of fracture for elliptic problems with flux and solution jumps[END_REF][START_REF] Ramière | A general fictitious domain method with immersed jumps and non-conforming structured mesh[END_REF]. Among other ones, well known methods of the second class of methods are the Immersed Boundary (IBM) [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF] and the Ghost Cell Methods [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF][START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF], the Fictitious Domain Method based on the explicit use of Lagrange multipliers [START_REF] Glowinski | A fictitious domain method for dirichlet problem and applications[END_REF], the Penalty methods [START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Ramière | A fictitious domain approach with spread interface for elliptic problems with general boundary conditions[END_REF][START_REF] Ramière | Convergence analysis of the Q1-finite element method for elliptic problems with non-boundary-fitted meshes[END_REF][START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF], the FAT boundary method [START_REF] Maury | A Fat Boundary Method for the Poisson problem in a domain with holes[END_REF] and the Finite Cell Method [START_REF] Duster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF]. These techniques greatly reduce the computational cost because re-meshing is in principle not required as there is only fluid cells. Concerning the IBM, this class of non-boundary conforming techniques has been first introduced by Peskin [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF] in order to simulate the incompressible blood flow in a human heart by adding a vectorial supplementary term, referred to as the forcing term. This continuous term is based on a Dirac delta function smeared over a stencil of few Cartesian nodes, in the continuous Navier-Stokes equations. Following Peskin's works about elastic structures, several IB-like methods for the Immersed Boundaries (IB), with different kinds of forcing terms have been proposed in the literature, such as the continuous Feedback Forcing (FF) of Goldstein & al. [19]. In this virtual-boundary method, the forcing term is spread over the immersed boundary and can be viewed as a force density. However, the FF term depends on the flow history and, due to large values of the stiffness coefficients, its application is limited when considering BC on rigid structures. The Direct Forcing (DF) method is an alternative approach proposed by Mohd-Yusof [START_REF] Mohd-Yussof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF] and then adapted by Fadlun & al. [21]. It consists in directly applying the desired boundary conditions on the Cartesian nodes in the fluid near the fluid-solid interface (through the characteristic function of the solid), leading to a sharp representation. Thanks to this formulation, the forcing term can be easily computed as it does not depend on the flow history [START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF]. The accuracy of this method lies in the direct forcing numerical scheme itself because the calculation of the forcing term is based on the discretized form of the governing equations. Since Mohd-Yusof, the DF method has been successfully applied to various FSI problems, cf. [START_REF] Kajishima | Interaction between particle clusters and particle-induced turbulence[END_REF]- [START_REF] Chiu | A differentially interpolated direct forcing immersed boundary method for predicting incompressible Navier-stokes equations in time-varying complex geometries[END_REF] and [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF]- [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF] for instance. As we would like to simulate complex FSI problems, such as the fluid-elastic instability of the heat-exchanger tube bundles, we have adopted the DF approach due to its simplicity of implementation and the low associated computational cost.

One undesirable property of IBM is the generation of spurious force oscillations (SFOs) [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF][START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF][START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF][START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF][START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF][START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] when dealing with moving bodies on a fixed computational grid. It is observed for all type of IBM, including discrete [START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF] and distributed [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF] forcing term methods or ghost-cell methods [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF][START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. SFOs degrade the quality of solutions and the computation of the stress forces near the fluid-solid interface. Therefore, whatever the considered IBM, one has to deal with these SFOs in order to simulate moving bodies under prescribed or flow induced movement. In the literature, it is reported that the main source of the SFO is the temporal discontinuity in the velocity at the grid points, called dead cells, where fluid becomes solid due to the body motion [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF][START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF]. The magnitude of the spurious force oscillations, mainly attributed to the pressure (through the velocity divergence), are said to decrease with the decrease of the space-step ∆x and the increase of the time-step ∆t as O( ∆x α ∆t ) with α a power-law coefficient [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF][START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. Magnitude of the SFO also decreases when using a carefully designed discrete delta function [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF]. Forcing into the solid (including dead cells) [START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF] are known to reduce the SFO without suppressing them. Few techniques have been proposed to reduce these spurious force contributions and can be classified in two categories. On the one hand, one preferentially works on the immersed boundary conditions of the momentum balance equation. With the solid volume fraction as weight, Kajishima & al. [START_REF] Kajishima | Interaction between particle clusters and particle-induced turbulence[END_REF] avoid the spurious pressure oscillations. Using a referential linked to the solid, Kim & Choi [START_REF] Kim | Immersed boundary method for flow around an arbitrarily moving body[END_REF] eradicate the SFO, but their technique is devoted to a unique solid. Luo & al. [30,[START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF] propose an implicit hybridization (through an iterative method) of the flow-solver numerical scheme and of the IBM in order to deal with the free-dead cell transitions and Chiu & al. [START_REF] Chiu | A differentially interpolated direct forcing immersed boundary method for predicting incompressible Navier-stokes equations in time-varying complex geometries[END_REF] replace the algebraically-interpolated method by a differentially interpolated method. On the other hand, one may prefer working on the immersed boundary conditions of the mass balance equation. Kim & al. [START_REF] Kim | An immersed-boundary finite volume method for simulations of flow in complex geometries[END_REF] and Lee & al. [34] combine a DF approach and an extra mass source term in the mass balance equation. This minimizes the pressure oscillations coming from the spatial discontinuities of the pressure near the IB when fresh cells (cells in the fluid and previously in the obstacle) are released. Seo & al. [36] conjointly use a DF ghost-cell method [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] for the Navier-Stokes equations and a cut-cell approach to design the pressure solver. Lee & You [START_REF] Lee | An implicit ghost-cell immersed boundary method for simulations of moving body problems with control of spurious force oscillations[END_REF] combine the same ghost-cell method and the Kim & al. approach. This work is twofold. First, it is devoted to the analysis of the SFO. Through numerical experiments and theoretical analysis, we confirm previous findings [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF][START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF] that designate the dead-cell contribution as the main source of the SFO. This is stated too from the confrontation of the SFOs and the fresh-cell apparition occurrences. We propose in this paper a more formal derivation of the magnitude order of the SFO as O( ∆x 2 ∆t ) + O(∆x) for a first-order in space interpolation. Using rigorous numerical evaluations based on Taylor expansions, we produce the specific contributions depending on the kind of considered cells (dead, fresh or unchanged). Secondly, this work provides a practical remedy to SFOs using a regularization technique for the DF method. In the philosophy of the hybridization approach promoted by Luo & al. [30,[START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF] or of the differentially interpolated method [START_REF] Chiu | A differentially interpolated direct forcing immersed boundary method for predicting incompressible Navier-stokes equations in time-varying complex geometries[END_REF], we realize an explicit weighted mixing between the interpolation scheme and a blind version (without obstacles) of the flow solver. In fact, our method differs from that of Luo & al. by an explicit forcing term (no iteration), leaving the native numerical scheme unchanged, and by the determination of the weight. The former property is very useful to introduce the DF method in previously existing industrial codes. Formally, it leads to a weighting of the forcing term that can be understood as the definition of a regularized characteristic function of the solid. Unlike [START_REF] Kajishima | Interaction between particle clusters and particle-induced turbulence[END_REF], this weighting is not restricted to the solid volume fraction only. Its design is mainly led by the search of a smooth transition between the Navier-Stokes velocities and the forced velocities in the vicinity of the IB. As in [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF], the method retains the sharp-interface representation of the solid body surface. In practice, we only modify the characteristic functions of the obstacles and consider the standard explicit DF method. Under CFL restriction, only one strip of fresh or dead cells can appear during a time step. All the cells inside the obstacles are forced. Moreover using a non-incremental projection scheme [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires[END_REF], we compute the pressure and the velocity correction through all the fictitious domain (fluid and solid). For clarity purpose, our presentation is done in 2-dimension space. But, that does not constitute a restriction of the method which is actually carried out in 3-D.

We first present in Section 2 the fluid governing equations, the numerical scheme, the DF method and, following Introïni & al. [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF], the interpolation schemes used to implement it with order one or two in space. Then, the SFO origins are numerically analyzed on the Seo and Mittal test case [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] in Section 3. The O( ∆x 2 ∆t ) + O(∆x) dependency, coming from the dead cells, is theoretically demonstrated in an original way for the first-order in space version of the standard DF. Following the results of this analysis, we propose in Section 4 what we call the regularized (RG) approach. This approach, easy to implement, significantly cuts off the SFOs without increasing the computational cost and conserves the accuracy and the space convergence order. With the same assumptions as before, we theoretically show that the magnitude order of the SFO varies as O(∆x + ∆t) in the case of the RG DF. Confirming our analysis, numerical and physical validations of the standard DF and RG DF approaches are conducted in Section 5 on various test cases ranging from the Taylor-Couette problem to a cylinder with an imposed sinusoidal motion subjected to a cross-flow.

Governing equations and numerical method

Governing equations

This paper deals with incompressible flows around moving solid obstacles. The full domain is named Ω = Ω F ∪ Ω S , where Ω F is the fluid domain and Ω S the solid one, and the fluid-solid interface is indicated by Σ, see Fig. 1.

The governing equations used to describe unsteady incompressible flows are given by: The knowledge of the local fluid force ρν(∇u + ∇ t u) • npn or of its integration Ω S (ρν(∇u + ∇ t u) • npn)dΣ on the fluid-solid interface Σ allows to access the dynamic of the body if not prescribed. In turn, this body motion induces a change in the flow structure. Although crucial for the fluid-structure interaction, this point is not discuss in this paper.

∂u ∂t + ∇ • (u ⊗ u) + 1 ρ ∇p -ν∇ 2 u = f in Ω F ∇ • u = 0 in Ω F (1

Spatial discretization and numerical method

The equations are discretized on a uniform Cartesian mesh using a finite volume approach with a staggered grid arrangement of the variables (u, P). As a result, the pressure degrees of freedom are located at the cell centers whereas those of each velocity component are placed at the middle of the cell edges as presented in Fig. 2. The governing Eqs. (1) are integrated over each control volume ensuring conservation of mass and momentum balance. Here, the convection and diffusion terms are respectively approached by the QUICK and centered schemes [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. A non-incremental fractional-step scheme is used to solve the incompressible Navier-Stokes equations. The fractional step method or projection method was introduced by Chorin and Temam in 1968-1969 for incompressible flows [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires[END_REF]. On the basis of their work, several variants have been proposed and the reader is referred to Guermond [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for a recent overview of these methods. In L ∞ norm, the rate of convergence in time of the non-incremental version is theoretically first-order [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF].

Direct Forcing Approach

As previously discussed in the introduction, there are several types of immersed boundary method (IBM) [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF][START_REF] Goldstein | Modeling Modeling a no-slip flow boundary with an external force field[END_REF][START_REF] Mohd-Yussof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF]. They all consist in computing flows around complex geometrical shapes (static/moving/deforming bodies) by immersing the physical domain Ω F in a simpler fictitious one Ω. The presence of embedded time-varying geometry Ω S in the computational domain Ω is taken into account thanks to a source term F IBM that is added in the Navier-Stokes equations [START_REF] Farhat | Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution[END_REF]. In this work, the domain Ω S is discretized through a 3-D Lagrangian triangulation. At the end of each time step, the Lagrangian mesh is moved following the fluid-solid interface velocity. Following the direct forcing (DF) method developed in [START_REF] Mohd-Yussof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF], a discrete source term F DF imposes immersed boundary conditions for the flow in the vicinity of the obstacles and is added in the discretized Navier-Stokes equations (here, fully explicit for the sake of the presentation) as follows:

u n+1 -u n ∆t + ∇ h • (u n ⊗ u n ) + 1 ρ ∇ h p n+1 -ν∇ 2 h u n = F DF in Ω (2)
∇ h • u n+1 = 0 in Ω (3)
where the subscript • h stands for discrete operators (omitted hereafter). The source term F DF refers to a direct forcing term [START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF] and is defined by:

F DF = χ m u m -u ∆t in Ω, ( 4
)
where u is the predicted velocity of the fractional-step scheme without IB, χ m is the characteristic function of the imposed velocities domain Ω m I , which is different depending on the considered model, and u m = u m (u S , u ) is the imposed velocity resulting from the adopted reconstruction model. Following the authors, the reconstruction model may be Lagrange interpolations following the grid directions [START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF] or in a multi-direction approach [START_REF] Kim | An immersed-boundary finite volume method for simulations of flow in complex geometries[END_REF], the normal direction to the body surface [START_REF] Yang | An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries[END_REF], power law interpolations [START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF] or least squares interpolations [START_REF] Peller | High-order stable interpolations for immersed boundary methods[END_REF]. In the rest of the paper, the subscript m means "model" that can be "base" (first order in space model) or "linear" (second order in space model).

The resolution of Eq. (2-3) is split according to the following algorithm:

Step 1: Computation of the predicted velocity u without taking into account the immersed boundary conditions:

u -u n ∆t + ∇ • (u n ⊗ u n ) -ν∇ 2 u n = 0 in Ω (5)
Step 2: Computation of the imposed velocity and the forcing term (see Sec. 2.3.1 and 2.3.2):

u m = u m (u S , u ) ; F DF = χ m u m -u ∆t (6)
Step 3: Addition of the direct forcing term:

ũ = u + ∆tF DF (7)
Step 4: Projection:

1 ρ ∆p n+1 = 1 ∆t ∇ • ũ in Ω (8)
Step 5: Correction:

u n+1 = ũ - ∆t ρ ∇p n+1 in Ω (9)
In the following section we detail the interface models used to calculate the imposed velocities.

The base model

The simplest approach to implement this IBM is to consider a first order accurate in space method. As for the ghost cell approach [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF], we choose to only force velocities that are inside the solid domain, depicted by red points in the Fig. 3. Due to this formulation, Ω base I = Ω S and χ base is defined as follows:

χ base (x) = 1 if d(x) ≥ 0 0 otherwise.
where d(x) is the signed distance function to the fluid-solid interface with d(x) ≥ 0 inside Ω S . At the beginning of the time step, the signed distance is computed in 3-D using the barycenters and the external normal vectors of the Lagrangian-mesh faces. In this formulation, the imposed velocity is the solid velocity:

u base (x) := u S (x). ( 10
)
With a Cartesian grid, the base model leads to a stepwise description of the immersed interface Σ.

The linear model

A way to improve the spatial order of the method is to use an interpolation scheme to reconstruct the velocity field near the interface. This is the most widely employed approach in the literature [START_REF] Fadlun | Combined immersed-boundary finite-difference for three-dimensional complex flow simulations[END_REF] and we have adopted the procedure developed by Introïni & al. [39]. The linear model consists in using an interpolation scheme to calculate the fluid velocities located near Σ. More precisely, Ω linear I is divided into two sub-domains as depicted in Fig. 3. The first one is made up of the velocities located inside Ω S that are forced to u S as for the base model. The second one must contain at least the fluid velocities having a neighbour in Ω S and are forced using an interpolation scheme to guaranty the second order global accuracy. Such a domain can be defined by adding a ∆x-thickness layer to the solid domain, leading to the following expression of χ linear :

χ linear (x) = 1 if d(x) ≥ -∆x 0 otherwise.
where ∆x is the cell length. Then, u linear takes the following form:

u linear (x) := u int (x) if -∆x ≤ d(x) < 0 u S (x) if d(x) ≥ 0 (11)
where u int is the velocity obtained with a linear interpolation scheme, second order accurate in space. In this interpolation scheme, the fluid contribution is built through an averaged reconstruction of the velocity gradient near Σ and the solid contribution is determined by means of a minimization problem, see [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF]. This model is illustrated in Fig. 3, where the red velocities are forced to the solid velocity (as for the base model) and the green one to interpolated velocities (these velocities are free when considering the base model).

Σ linear I = {x, d(x) = -∆x} Σ = Σ base I = {x, d(x) = 0}

An analysis of the Spurious Force Oscillations

This part is dedicated to the analysis of the origin of the SFO and to the mathematical evaluation of its magnitude. First of all, we introduce the SFO through the Seo and Mittal test case [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. As already mentioned, without restriction on the space dimension d, dependency on O( ∆x 2 ∆t ) was notified by Lee & al. [34] through a second-order interpolation. Considering only the mass balance equation, Seo and Mittal theoretically demonstrated dependency on O( ∆x d ∆t ) [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. In this work, we mainly focus on the momentum balance equation. Then using Taylor expansions of the forcing term and of the velocity near the IB, we quantify each specific contribution to the SFO magnitude coming from the dead, fresh or unchanged cells. Finally numerical studies give confirmation of this power-law.

Seo and Mittal test case

In order to introduce the issue of the generation of spurious force oscillations, when using a DF approach for a rigid moving body, the simple test case of Seo and Mittal [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] is considered. It consists in a circular cylinder centered at the coordinates (x c , y c ) of diameter D which oscillates sinusoidally along the x-direction in a fluid at rest :

x c (t) = x c (0) + X 0 (1 -cos(2π f 0 t)) y c (t) = y c (0)
where (x c (0), y c (0)) are the initial coordinates of the cylinder center, X 0 is the amplitude of the oscillation and f 0 is its frequency (see the Fig. 4(a) for a schematic representation). The period of oscillation is T 0 = 1 f 0 . The computational domain is a 4D × 4D square with the upper and lower boundaries set to no-slip walls, homogeneous Neumann conditions for the pressure and the velocity on the left and right sides, (x c (0), y c (0)) = (2D, 2D) and X 0 = 0.125D. Due to the cylinder motion, the cylinder maximum velocity is U 0 = X 0 2π f 0 , the Reynolds and Strouhal numbers, respectively defined by Re = U 0 D ν and S t = f 0 D U 0 are set to 78.5 and 1.27. In SFO studies, we are interested in the temporal variations of the following physical dimensionless coefficients: the pressure drag C P , the friction drag C ν and the total drag C T . They are the components along the x-axis of the pressure, viscosity and the total force. In addition, we consider the total contribution of the direct forcing term C DF .

These coefficients are defined as follow:

C P = Ω S -pndΣ 0.5ρD 3 f 2 0 (12) C ν = Ω S ρν(∇u + ∇ t u) • ndΣ 0.5ρD 3 f 2 0 (13) C T = C P + C ν (14)
C DF = Ω S ρF DF • e x dΩ 0.5ρD 3 f 2 0 ( 15
)
By integrating the Navier-Stokes equations including the forcing term, C DF is linked to C T by the following expression:

C T = -C DF + Ω S ρ ∂u ∂t + ∇ • (u ⊗ u) dΩ • e x 0.5ρD 3 f 2 0 . ( 16
)
The volume space integration of Eq. ( 15) or Eq. ( 16) is numerically done on the imposed-velocity domain Ω m I and not on the exact solid domain Ω s . We use the velocity control volumes of our FV scheme, cf. Fig. 2. The surface space integration of Eq. ( 12) or Eq. ( 13) is not done on the exact interface ∂Ω s but on the approximated fluid-solid interface Σ m I using the velocity-control-volume faces included in Σ m I . As the pressure and the Cauchy stress are defined on these faces, we directly use the terms involved in the momentum balance equation.

To illustrate the SFO phenomena, the problem is first solved on a 64 2 uniform grid with a fixed time-step ∆t = 0.002 1 f 0 = 0.002T 0 , corresponding to the following CFL number 0.025, where CFL = U 0 ∆t ∆x and ∆x is the grid spacing. Considering this grid, the cylinder diameter is spanned by 16 grid points. The computation is performed with the DF method described in subsection 2.3. The evolution over a period T 0 of the dimensionless coefficients of interest are presented in Fig. 4(b). It is important to point out that it is mainly the pressure drag component that is disturbed by the spurious oscillations whereas the time history of the friction drag is almost regular. Therefore, the SFOs are in fact due to the spurious pressure oscillations (SPO). Then, it is also interesting to notice that C DF is nearly equal to C T .

Then, the convergence in space is investigated by considering a fixed time-step ∆t = ∆t 0 = 0.002T 0 , with four uniform grids (64 2 , 128 2 , 192 2 and 256 2 cells), leading to the following CFL numbers: 0.025, 0.05, 0.1 and 0.2. Thus, on these four grids, the cylinder diameter D is spanned by 16, 32, 48 and 64 grid points, respectively. The convergence in time is investigated by considering a 64 2 -cell grid and the following time-steps ∆t = ∆t 0 , 2∆t 0 , 4∆t 0 and 8∆t 0 , leading to the following CFL numbers: 0.025, 0.05, 0.1 and 0.2. Figs 5(a) and 5(b) display the C P time histories over a period of time. Qualitatively, they demonstrate that the pressure oscillations decreases with decreasing grid spacing and with increasing time-step, as mentioned in [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF][START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF].

Theoretical analysis of the SFO origin

We analyze the SFO origin in the case of the DF base model described in subsection 2.3.1. First of all, we assume the origin is located in the forcing step (2. and 3.) between the prediction (1.) and the projection steps (4.) of this algorithm. Indeed, due to the nature itself of the DF approach, the projections (4) and correction steps (5.) are blind to the solid boundary as they do not take directly into account the solid boundaries for their computation. However, the projection step spreads the pressure oscillations over all the computational domain via the source term of the Poisson equation which can be disturbed because of the forcing step.

As previously shown on Seo and Mittal's test case, C DF and C P are disturbed by equivalent spurious oscillations. In order to quantify the amplitude of these oscillations, it is assumed a constant solid velocity in space and time (the analysis remains valid for regular variations). This assumption allows Taylor expansions on the velocity field in time and space. For the sake of simplicity, we assume a CFL number smaller than unity and the following hypothesis: where Ω n s,h is the set of velocities forced at the time-step n. Eq. ( 17) means that all the velocities forced at the time-step n were already forced at the time-step n-1 even if that usually does not hold for two consecutive steps in real life. But it enables a mere simplification of the analysis consistent with a limited-space presentation. However, it is not essential to the computations as it will be stated below. This reduces the inventory of the different cells categories following their history with respect to the forcing term. The velocity evolution between the time-steps n and n + 1 is schematically described as follow:

Ω n-1 s,h = Ω n s,h Ω n+1 s,h (17)
u n ∆tRHS n (u n ) --------→ O(∆t) u ,n ∆tF DF ----→ O(∆t) ũ -∆t ρ ∇p n+1 -------→ O(∆t) u n+1 ∆tRHS n+1 (u n+1 ) -----------→ O(∆t) u ,n+1 (18)
where ∆tRHS n (u n ) denotes the explicit Navier-Stokes solver of the prediction step (1.). One can measure the DF term variation, between the time-steps n and n + 1, by calculating:

Ω n+1 s,h F n+1 DF dΩ - Ω n s,h F n DF dΩ = Ω n+1 s,h ∩Ω n s,h F n+1 DF -F n DF dΩ+ (19) Ω n+1 s,h Ω n s,h F n+1 DF dΩ - Ω n s,h Ω n+1 s,h F n DF dΩ.
By expressing the integrals at the discrete level and writing the DF term expression, the above equation becomes:

Ω n+1 s,h F n+1 DF dΩ - Ω n s,h F n DF dΩ = K∈Ω n+1 s,h ∩Ω n s,h u S -u ,n+1 -u S + u ,n ∆t ∆x d + ( 20
)
K∈Ω n+1 s,h Ω n s,h u S -u ,n+1 ∆t ∆x d - K∈Ω n s,h Ω n+1 s,h u S -u ,n ∆t ∆x d
where K is a mesh element and ∆x d the volume of a cell in d dimension. Denoting FC (fresh cell) the set of velocities that became freshly fluid (K ∈ Ω n s,h Ω n+1 s,h ) and DC (dead cell) the set of velocities that became freshly forced (K ∈ Ω n+1 s,h Ω n s,h ), we get:

Ω n+1 s,h F n+1 DF dΩ - Ω n s,h F n DF dΩ = ∆x d ∆t            K∈Ω n+1 s,h ∩Ω n s,h u ,n -u ,n+1 + (21) DC u S -u ,n+1 - FC u S -u ,n        .
Forced Cells at n and n + 1. u ,nu ,n+1 is associated to the velocities forced at time-steps n and n + 1. Thanks to Eq. ( 18), u ,n+1 = u n+1 + O(∆t) = ũn + O(∆t) = u S + O(∆t). By doing the same analysis and thanks to the assumption [START_REF] Duster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF], u ,n = u S + O(∆t). Therefore:

u ,n -u ,n+1 = O(∆t). (22)
Dead Cells. u Su ,n+1 is associated to velocities free at n and forced at n + 1. As this velocity is forced at n + 1, it has at least one forced neighbor at the previous time-step n as we assumed CFL < 1. In addition, the DF base model is of order one in space, implying that u ,n+1 = u ,n+1 Forced neighbor + O(∆x) where u ,n+1 Forced neighbor = u S + O(∆t). Thus:

u ,n+1 -u S = O(∆x) + O(∆t). ( 23
)
Fresh Cells. u Su ,n is associated to velocities forced at n and free at n + 1. As these velocities were previously forced at n -1 due to the assumption 17, u ,n = u n + O(∆t) = u S + O(∆t). Thus:

u S -u ,n = O(∆t). ( 24
)
Conclusion. Finally, the number of FCs and DCs increases with the mesh size. These cells appear by a layer of thickness adjacent to the boundary (dimension 1 in 2D and dimension 2 in 3D). Their number is proportional to

( D ∆x ) d-1 = O( 1 ∆x d-1
). The number of cells forced at the time-steps n and n + 1 increases at the same rate as the mesh elements proportional to

( D ∆x ) d = O( 1 ∆x d ).
As last with the assumption 17:

Ω n+1 s,h F n+1 DF dΩ - Ω n s,h F n DF dΩ = ∆x d ∆t O( ∆t ∆x d ) + O( ∆t + ∆x ∆x d-1 ) + O( ∆t ∆x d-1 ) = O(1) K∈Ω n+1 s,h ∩Ω n s,h + O(∆x) + O( ∆x 2 ∆t ) K∈DC + O(∆x) K∈FC (25) = O(1) + O(∆x) + O( ∆x 2 ∆t ).
This analysis can also be conducted without the assumption [START_REF] Duster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF]. But it would have been necessary to distinguish the velocities free at n -1 that are forced at n and the velocities already forced at n -1. However for the sake of clarity this development is not exposed as the conclusions are exactly the same.

As the value of the term O( ∆x 2 ∆t ) can be big, specially for small ∆t, we can guess that the DCs play the principal role in the SFO. But as the coefficients of the various contributions are unknown, we cannot confidently predict whether the DC or the FC are the main contribution to the SFO. To remove this uncertainty Fig. 6 displays the time history of the forcing term maximum, the number of DCs and FCs with the numerical conditions of Section 3.1. Fig. 6(a) shows that the forcing term oscillations occur at the same times as new cells enter in the solid domain and become dead cells. Whereas, on Fig. 6(b) it can be noted that the fresh cells do not trigger significant oscillations as no time coincidence can be found between the FC's occurrence and the forcing term maximum. This means the DC are the main source of pressure oscillations (at least for small time-steps). The FC contributions are commonly attributed to space discontinuities of the velocity and the pressure [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF]. Here the non-incremental pressure-projection scheme and the computation of velocities and pressures over the whole fictitious domain should contribute to reduce this space discontinuity for FCs.

Then thanks to Eq. ( 25), for small values of the time-step and the DF method with a first-order interpolation, SFOs are proportional to O( ∆x 2 ∆t ). Through a more formal derivation, giving the specific contributions of each kind of cells (dead, fresh or unchanged), we generalize here the Lee & al. expression O( ∆x 2 ∆t ) [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF], based on simple considerations concerning the velocity change for dead cells.

Theoretical analysis of the SFO magnitude using the linear DF was not conducted here. But the numerical experiments described in Section 4.3 conduct to a power-law about O( ∆x 3 ∆t ) for small time-steps.

Numerical validation on SFO

As above mentioned, Eq. ( 25) is consistent with the expression O( ∆x α ∆t ) found in the literature. However, it can be interesting to point out that the biggest is the time-step the more the oscillations tend to be dependent on O(∆x) instead of O( ∆x 2 ∆t ). Taking into account all the specific contributions, the present analysis brings an extra light into the behavior of the SFO amplitude. This point is numerically illustrated here. The pressure 2δ-discontinuity, noted C 2δ P and introduced by Seo and Mittal [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF], is used to obtain quantitative results. It is defined by the following equation:

C 2δ P = C n+1 P -2C n P + C n-1 P ( 26
)
where n is the time-step index.

The result obtained for C 2δ P,Max is in very good agreement with Eq. ( 25). In Fig. 7(b), it is interesting to notice that when the time-step increases (corresponding at the left side of the picture), C 2δ P,σ and C 2δ P,Max tend to decrease slower with the time-step than for the smallest time-step (right side of the picture). This original result is well described by our theoretical analysis summarized by Eq. ( 25), because the contribution of the DC is O(∆x) + O( ∆x 2 ∆t ). So when ∆t increases, O( ∆x 2 ∆t ) decreases and the contributions in O(∆x) of the DC and FC become significant in comparison with O( ∆x 2 ∆t ). At the opposite, when the time-step is very small the contribution in O( ∆x 2 ∆t ) of the DC is much higher than the other ones.

C 2δ P,σ is quite close to the theoretical result also. It is necessary to have in mind that the maximum of C 2δ P quantify mainly the oscillations triggered by the DC. This explains the behavior difference between C 2δ P,σ and C 2δ P,Max . Indeed, as the main source of oscillations, the DC are at the origin of the biggest peaks (which are the maximums), whereas the FC contribution is much lower. It is why C 2δ P,Max is rather close to the predicted DC variation O( ∆x 2 ∆t ). The standard deviation measures the dispersion of C 2δ P . Thus, C 2δ P,σ quantifies all the oscillations contributions (FC and DC). Therefore, it cannot have a variation following exactly the one of the maximums.

The direct forcing regularized formulation

Analysis of Section 3 done with the base model predicts a SFO behavior in O( ∆x 2 ∆t ) + O(∆x) due to dead cells. Considering coarse to moderate grids and small to moderate time-steps (precisely the configuration where the SFOs are annoying) the power-law is O( ∆x 2 ∆t ). Numerical studies give confirmations of this power-law. A simple way to tackle the oscillations would be to decrease the grid-spacing and increase the time-step. However, that would be at the expense of the computational time when decreasing the grid-spacing and at the expense of the precision when increasing the time step. Consequently it cannot be an end in itself. To reduce the DC contributions to the spurious pressure oscillations, the mass balance can be corrected (e.g. [START_REF] Kim | An immersed-boundary finite volume method for simulations of flow in complex geometries[END_REF][START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]) or the imposed velocities can be carefully computed taking into account the Navier-Stokes solver (e.g. [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF]). These methods (cut-cell method for mass balance, fully-implicit hybridization of interpolated and Navier-stokes velocities, . . . ) are efficient but quite invasive considering for instance an industrial computer code already using a standard DF method. That is why a regularized (RG) formulation of the direct forcing, called the RG DF approach, easy to implement and capable of cutting off efficiently the oscillations without any extra computational time, is presented. Through a regularization process, it is possible to reduce the SFO by roughly an order of magnitude in comparison with the original DF approach. Moreover, the O(∆x) + O( ∆x 2 ∆t ) dependency of the SFO amplitude, with first-order in space interpolation, is changed to O (∆x) + O (∆t). This change in the power-law is very interesting for small timesteps. To the knowledge of the authors, it was not achieved by the past (e.g. the power-law is not changed when using the method described in [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]). First, we focus on the numerical method and analyze its impact on the spurious force oscillations. Then, we numerically study the space and time convergence of the oscillation amplitude on the Seo and Mittal's test case using the RG DF.

Numerical Method

In the standard DF method, the added force F DF is defined through the use of a Heaviside-like generalized function for χ m , cf. Eq (4), resulting in a "all or nothing" method and a sharp transition of the forcing term through the free-to-forced velocity interface. Moreover, when a free-velocity grid point is crossed by this interface, unless a specific treatment, no consistency in time is guaranteed between the previous Navier-Stokes solver value and the new interpolated one, contributing to the SFO apparition. Luo & al. [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF] and Chiu & al. [32] propose such kind of treatment to improve this time consistence for a ghost cell method, taking into account the Navier-Stokes solver in the design of the interpolation scheme. Here, the key idea of the RG DF formulation is, for a new given forced-velocity point, to smooth the transition between the values given by the Navier-Stokes solver and by the interpolation (or affectation) scheme.

Considering a fluid forced-velocity point in Ω m I , we propose the following new formulation for the discretized imposed velocity, now denoted by u i to distinguish from the interpolated one u m :

u i = ξ m u m + (1 -ξ m )u F RG = u i -u ∆t (28)
where ξ m (x) is a regular function of the signed distance d(x). This formulation is named "regularized" (RG) because, for the degrees of freedom with 0 < ξ m < 1, the new imposed velocity u i is regularized through a linear combination of the blind Navier-Stokes solver velocity u and of the interpolated velocity u m . We choose to define ξ m (x) in such a way that it varies linearly from 0 to 1 across a layer of thickness ∆d centered on the iso-distance d m 0 . But different functions ξ m (x) can be addressed. In this work, d m 0 is defined as the signed distance of the forced-free interface Σ m I (depending of the considered model m) and ∆d = ∆x. Hence the function expression can be interpreted as the forced-velocity domain fraction τ m (x) in the control volume of the forced-velocity point and reads as

ξ m (x) = ∆x + 2[d(x) -d m 0 ] 2∆x ≡ τ m(x) . (29)
This is graphically represented on Fig. 8 ≡ τ(x). This phase-mean approach can be found in [START_REF] Kajishima | Interaction between particle clusters and particle-induced turbulence[END_REF] for a first-order in space direct forcing method. For the linear model, the signed distance of the forced-free interface Σ linear I is d linear 0 = -∆x, cf Fig. 3. Then, the function ξ linear is similar to the function ξ base shifted by -∆x.

(x) = ∆x+2[d(x)-d m 0 (x)] 2∆x ≡ τ m (x) in 1D.
From a continuous-geometry point of view, there is no difference between the standard DF method and the regularized one.From a discrete-geometry point of view, the difference lies in the discretization of the indicator of the imposed-velocity region, leading to the corresponding discretized forcing term F RG :

F RG = ξ m u m -u ∆t . ( 30
)
that is the standard DF forcing term F DF with the generalized function χ m replaced by the regular function ξ m . In particular, F RG has exactly the same expression as F DF of the original formulation when the control volume of the fluid velocity is entirely in the forced domain (i.e.: ξ m = χ m = 1). In a certain way, considering this ξ m -weighting of the force, we can found some links with the volume-weighting of Kajishima and Takiguchi [START_REF] Kajishima | Interaction between particle clusters and particle-induced turbulence[END_REF] or the Uhlmann method [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] leading also to a weighting of the discrete force term on the Cartesian grid. But, in the proposed regularized DF, the weighting is not restricted to the only solid volume fraction and a sharp interface is considered. In another way, similarities can be found with the Luo & al. method [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF] in which an implicit hybridization of the flow-solver numerical scheme and of the velocity interpolation scheme is done. But here, this hybridization is explicit in time and it saves us to modify the discrete operators to include the interpolated velocity.

Practically, the implementation of Eq. ( 28) is:

u i =          u m if ξ m = 1 ξ m u m + (1 -ξ m )u if 0 < ξ m < 1 u if ξ m = 0. ( 31
)
Thereby, when the solid is in motion, the fluid velocity will be forced progressively (and not abruptly as in the original DF approach) to the interpolated forced velocity u m as it enters in the forced domain Ω m I , smoothing the temporal discontinuity of the velocity through the interface. This procedure is illustrated on Fig. 9(a) with ξ m = τ m . The black dashed-line square denotes a control volume occupied by a fraction τ n m (in grey) of the forced domain. Thus, the correspondent imposed velocity is proportional to this fraction. At the next time-step, cf Fig. 9(b), the forced domain has been vertically translated and the control volume is now occupied by a more consequent fraction τ n+1 m > τ n m . Thus, the imposed velocity is forced proportionally to this new forced domain fraction. All the other steps of the resolution of the incompressible Navier-Stokes equations are solved as previously exposed. As anyone can notice, the above formulation does not increase the computational time cost and is very easy to implement.

Finally, let us focus on the cut-cell mass/volume conservation with the RG approach. Including the definition (28) of the imposed-velocity u i in the discrete ũ-velocity divergence of the cut cell K,

∇ • ũ| K = 1 |K| f | f | u i f • n f , we get: ∇ • ũ| K = 1 |K| f /ξ m f =0 | f | u f • n f + 1 |K| f /ξ m f >0 | f | (ξ m f u m f + (1 -ξ m f ) u f ) • n f (32)
where |K| stands for the measure of the cut cell and f for a face of the cell of measure | f |. Then, this expression can be recast into a no-weighting contribution and a correction q:

∇ • ũ| K = 1 |K| f /ξ m f =0 | f | u f • n f + f /ξ m f >0 | f | u m f • n f -q (33) q = 1 |K| f /ξ m f >0 | f | (1 -ξ m f ) (u m f -u f ) • n f . ( 34
)
The term with the brackets in the RHS of (33) corresponds to ∇• ũ| K in case of the standard DF. Hence, we can consider that the RG formulation formally modifies the RHS of the pressure equation ( 8) by adding the source/sink term q. It differs from the Kim & al. correction [START_REF] Kim | An immersed-boundary finite volume method for simulations of flow in complex geometries[END_REF] on the solid velocity. Considering ξ m = τ m , (34) can be viewed as a fluid correction proportional to the fluid fraction (1 -τ m f ) of the cut-cell faces.

Impact on the oscillations

In order to evaluate the temporal variation of the forcing term with the RG formulation, we repeat the analysis conducted on the base model in the subsection 3.2 by replacing F DF by F RG in Eq (19)

Ω n+1 s,h F n+1 RG dΩ - Ω n s,h F n RG dΩ = Ω n+1 s,h ∩Ω n s,h F n+1 RG -F n RG dΩ+ (35) Ω n+1 s,h Ω n s,h F n+1 RG dΩ - Ω n s,h Ω n+1 s,h F n RG dΩ
taking into account ξ base ≡ τ in Eq [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF]. Ω n s,h is the set of velocities forced at the time step n. It can be divided into two subsets. The first one is the subset of the velocities with τ n = 1 and the second one is the subset of the velocities with 0 < τ n < 1. Again Ω n+1 s,h ∩ Ω n s,h is the set of velocities forced at the time-steps n and n + 1 and Ω n+1 s,h

Ω n s,h is the set of velocities that were free at the time-step n (τ n = 0) but that are forced at the time-step n + 1 (0 < τ n+1 ≤ 1). Similarly Ω n s,h Ω n+1 s,h is the set of velocities that were forced at the time-step n (0 < τ n ≤ 1) and that become free at the time-step n + 1 (τ n+1 = 0).

By expressing the integrals at the discrete level and by writing the regularized forcing term expression, one gets:

Ω n+1 s,h F n+1 RG dΩ - Ω n s,h F n RG dΩ = K∈Ω n+1 s,h ∩Ω n s,h τ n+1 u S -u ,n+1 -τ n u S -u ,n ∆t ∆x d + ( 36
) K∈Ω n+1 s,h Ω n s,h τ n+1 u S -u ,n+1 ∆t ∆x d - K∈Ω n s,h Ω n+1 s,h τ n u S -u ,n ∆t ∆x d
Using the notations introduced before, it comes:

Ω n+1 s,h F n+1 RG dΩ - Ω n s,h F n RG dΩ = ∆x d ∆t            K∈Ω n+1 s,h ∩Ω n s,h τ n+1 u S -u ,n+1 -τ n u S -u ,n + (37) DC τ n+1 u S -u ,n+1 - FC τ n u S -u ,n        .
Dead Cells. τ n+1 u Su ,n+1 is associated to the dead cells. That means that the velocity was free at the time-step n (τ n = 0) and is forced for the first time at the time-step n + 1. We consider the most penalizing case where the fraction of forced domain in the control volume is maximum while it was zero at the previous time step. This occurs when d = -∆x/2 at time n, see Fig. 10(a). Indeed, the distance at time n + 1 is d = -∆x/2 + |u S |∆t, see Fig. 10(b). Considering τ = ∆x+2d 2∆x (see Fig. 8), we obtain τ n+1 = U s ∆t ∆x ∝ CFL. A freshly forced velocity has at least one forced neighbor (located at a few ∆x) entirely forced at the previous time-step. For the base model (first order in space), this implies that: u ,n+1 = u ,n+1 Entirely f orced neighbor u S +O(∆t) +O (∆x). Then:

τ n+1 u S -u ,n+1 = O ∆t ∆x (O (∆x) + O (∆t)) .
Fresh Cells. τ n u Su ,n is associated to the fresh cells. That means the velocity was last forced at time-step n and will be free at n + 1. There are three possible histories for the FC. Either there exists a time-step n -N where this velocity was entirely forced (τ n-N = 1). Or this time-step does not exist and there remains two possibilities: either it has always been partially forced, either there exists a time-step n -N where it has been a DC.

Let's consider the first possible history. That means u ,n-N+1 = u S + O (∆t) by using Eq [START_REF] Maury | A Fat Boundary Method for the Poisson problem in a domain with holes[END_REF]. We reconstruct the history as follows:

                       u ,n -u S = 1 -τ n-1 u ,n-1 -u S + O (∆t) ... = ... u ,n+2-N -u S = 1 -τ n-N+1             u ,n-N+1 -u S O(∆t)             + O (∆t)
and we get: As τ < 1 due to its definition, N-1 i=1 1 -τ n-i < 1 = O (1) and 1

u ,n -u S = O (∆t) N-1 i=1 1 -τ n-i + O (∆t)         1 + N-3 k=0 k j=0 1 -τ n-1-j         .
+ N-3 k=0 k j=0 1 -τ n-1-j < (1 + (N -2)) = O (1), then: u ,n -u S = O (∆t) .
Now, let's consider the second and third possibilities. As this velocity has been forced at n, it has at least one neighbor (at a distance of a few ∆x) entirely forced at n. Thanks to the fact the base model with the RG approach is first order in space, we have:

u n = O (∆x) + u n Entirely f orced neighbor u S +O(∆t)
and therefore by using Eq [START_REF] Maury | A Fat Boundary Method for the Poisson problem in a domain with holes[END_REF]:

u ,n -u S = O (∆x) + O (∆t) .
Finally, we conclude that u ,nu S = O (∆t) + O (∆x) for every FC. Regarding τ n and using again the CFL condition, we found that τ n is O U s ∆t ∆x by similar considerations to those made for the Fig. 10 ( U s ∆t∆x ∆x 2 is the maximum added or removed solid fraction of a control volume in one time-step). Eventually, we get:

τ n u ,n -u S = O ∆t ∆x (O (∆t) + O (∆x)) .
Forced Cells at n and n + 1. τ n+1 u Su ,n+1 -τ n u Su ,n is associated to the velocities forced at the time-steps n and n + 1. The set Ω n+1 s,h ∩ Ω n s,h can be divided into three distinct subsets. The first one is the subset of the velocities with τ n-1 = τ n = τ n+1 = 1. The second subset is composed of the velocities with 0 < τ n-1 < 1, 0 < τ n ≤ 1 and 0 < τ n+1 ≤ 1. Finally the third subset is composed of the velocities that were fluid at n -1 and freshly forced at n (that are in fact the DC at n).

For the first subset, we are exactly in the same case as for the DF base model: u Su ,n+1u Su ,n = u ,nu ,n+1 . Consequently, the analysis previously conducted is still valid: u ,nu ,n+1 = O (∆t).

Let's consider the second subset. We apply exactly the same analysis as for the FC: u ,nu S = O (∆t) + O (∆x) and then we get:

τ n+1 u S -u ,n+1 -τ n u S -u ,n = O ∆t ∆x (O (∆t) + O (∆x)) .
Finally the velocities of the third subset are DC at the time-step n. As previously done, we get: u ,nu S = (O (∆t) + O (∆x)) that leads to:

τ n+1 u S -u ,n+1 -τ n u S -u ,n = O ∆t ∆x (O (∆t) + O (∆x)) .

Conclusion.

For the same reasons previously exposed, concerning the number of DF and DC, we obtain:

Ω n+1 s,h F n+1 RG dΩ - Ω n s,h F n RG dΩ = ∆x d ∆t O ∆t ∆x d + O ∆t ∆x (∆t + ∆x) ∆x d-1 + O ∆t ∆x (∆t + ∆x) ∆x d-1 + O ∆t ∆x ∆t ∆x d-1 = O (1) + O (∆t) + O (∆x) K∈Ω n+1 s,h ∩Ω n s,h + O (∆t) + O (∆x) K∈DC + O (∆t) K∈FC = O (1) + O (∆t) + O (∆x) . ( 38
)
Comparing this relation to Eq 25, one can notice that the term associated to the DC no longer depends on 1 ∆t and that the dependency on ∆x is linear instead of quadratic.

Numerical validation on SFO

As in the previous subsection 3.3, the impact of the RG formulation is investigated on the convergence in space and time of the SFO with the Seo & Mittal's test case [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. Here we consider the base and linear models.

Regularized base model

The results of the RG base model in comparison of the standard base model for the DF are displayed on Fig. 11. It is quite obvious to notice the dramatic decrease of the SFO by using the RG formulation instead of the original approach.

As previously, C 2δ P is used to quantify the SPO amplitude, cf. Eq [START_REF] Yang | An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries[END_REF]. Fig. 12 shows |C 2δ P | plotted as a function of grid spacing and time-step for both the RG and standard DF base models. It is noted that there is at least a difference of one order of magnitude between both the formulations at a given time-step and grid spacing. Then, C 2δ P,σ and C 2δ P,Max are plotted as a function of grid spacing in Fig. 13 and time-step in Fig. 14. We denote these quantities by C p 2δ ,σ and C p 2δ

,Max where stands for the standard or regularized version. When considering time and space variations around ∆t = ∆t 0 , cf. Fig. 13(a), and D/∆x = 16, cf. Fig. 14(a), they can be fitted by the following power-law:

C p 2δ

RG,σ = ∆t 0.377 ∆x 0.9493 and C p 2δ RG,Max = ∆t -0.0055 ∆x 1.1782 .

In subsection 3.3, we concluded that, in the configurations where the SFO can have a large amplitude (coarse grid spacing, little time-steps), the main term is of order O ∆x 2 ∆t and well represented by C p 2δ DF,Max . Here, we expect that this dependence will be changed in O ∆t ∆x ∆x 2 ∆t = O (∆x) due to the introduction of τ ∝ ∆t ∆x . This dependence is confirmed concerning the space convergence on Figs 13(a) to 13(c) for various fixed timesteps, and concerning the time-step convergence on Figs 14(a) and 14(b) for various fixed grid spacing. Indeed, the pressure peaks, measured by C 2δ P,Max , have no more dependence in the time-step when the grid spacing is coarse, explaining the null time-step and the one spatial order convergences for C p 2δ RG,Max . Concerning C p 2δ RG,σ , the spatial convergence order is conform to the analysis, but the time-step order do not follow particularly Eq. [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF]. However, Eq. [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF] shows a contribution proportional to O (∆t) leading to a time-step dependence of C p 2δ RG,σ as it measures all the contributions (FC and DC) to the SPO. This is the case in Fig. 14(a) -14(b), even if the order is lower than one precisely because all the contributions are not proportional to O (∆t).

Regularized linear model

This subsection is devoted to the study of the impact of the linear-interpolation scheme on the pressure spurious oscillations. Even tough this procedure enables to gain one order of numerical rate of convergence in space leading in theory to more precise physical results for a same mesh refinement, it is necessary to assess its impact on the SFO by using once again the Seo and Mittal test case.

The results are presented on the Fig. 15 for the standard and the RG versions of the DF. It can be noted that the linear interpolation procedure leads to C p 2δ DF,linear,Max values smaller than the C p 2δ DF,base,Max ones. As for the base model, the RG process decreases the C 2δ P,Max values for the linear interpolation also. But this regularization is far away more efficient for the base model than for the linear model. Hence, the C p 2δ RG,base,Max values are the smallest ones. We can notice that, for the linear interpolation scheme, we obtain space convergence orders of C 2δ P,Max higher than for the base scheme. This stands for the standard and the regularized DF. Without a theoretical analysis as done for the base scheme, we can only conjecture a link with the fact that the standard and RG DFs equipped with a linear formulation are second order accurate in space, cf. Section 5.1.

Δ t =0.00025T Δ t =0.0005T Δ t =0.001T Δ t =0.002T 1e
Δ t =0.00025T Δ t =0.0005T Δ t =0.001T Δ t =0.002T

Conclusion

By comparing our results to the literature and particularly to Fig. 14 in [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF], it is found that the C 2δ P,σ is decreased by roughly one order of magnitude when using the standard DF base model (Fig. 13(a) and 14(a)) in comparison with their original method. Then, using the RG base model enables to gain at least one (or two; depending on the space and time resolutions) additional order of magnitude leading to similar level of SFO in comparison with their cut cell method. The main difference lies in the implementation as the RG base model is far more easy to implement than a cut-cell method. Also, let us notice that the present results compare well with the ones of Lee & al. (Fig. 15 in [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF]). They consider only the space-step evolution of the SFO's RMS. For common ∆X/∆X 0 ∈ [0.16; 0.64], our RG base-model results (C 2δ P,σ , Fig. 13) as well as our RG linear-model results (C 2δ P,σ ≤ C 2δ P,Max , Fig. 15) are in the same range of values as that of Lee & al. (or lower).

Numerical and physical validation

This section is devoted to the numerical and physical validation of the RG base and linear models and the comparison with the results of the standard DF ones. Three cases have been selected by increasing gradually the complexity of the physical phenomena at stake from a laminar cylindrical Couette flow to a cylinder with an imposed sinusoidal motion subjected to a cross-flow.

Numerical Validation

Our goal is to assess the numerical rate of convergence of the standard and the proposed RG DF methods without (base model) and with the interpolation scheme, by doing a grid convergence study. A cylindrical Couette flow is considered and the results compared to the analytic solution by calculating the relative L 2 (Ω f ) norm 2 of the error given by:

2 = N f i=0 (u i -u i re f ) 2 N f i=0 u i re f 2 ( 40
)
where the superscript i denotes the i-th face of the Eulerian grid, u i re f a reference velocity calculated on the i-th face thanks to the analytic solution and N f the total number of faces in the fluid region Ω f . This problem has already been considered in the frame of IBMs [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF][START_REF] Peller | High-order stable interpolations for immersed boundary methods[END_REF][START_REF] Cheny | The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties[END_REF]. It allows to numerically estimate the order of spatial accuracy of the proposed methods on problem involving rotating geometries.

Fig. 16(a) presents the geometrical features of the computational domain Ω = [0, L] × [0, L] composed of the solid and fluid domains: Ω = Ω s ∪ Ω f . The immersed boundaries Σ 1 and Σ 2 mimic the inner and outer cylinders of radius r 1 and r 2 . The inner cylinder rotates clockwise (ω 1 > 0) while the outer cylinder rotates in the counterclockwise (ω 2 < 0). In all calculations, we assume the Reynolds number Re = 40). A quasi-linear numerical rate of convergence is obtained for the base model and a quasi-quadratic rate for the linear interpolation scheme for both the DF and the RG formulations. These results confirm the enhancement of accuracy expected with the linear interpolation model, and show that the regularization method do not deteriorate the convergence order.

Physical Validation

A static cylinder in a cross-flow with Re = 100 and an imposed sinusoidal-motion cylinder subjected to a crossflow of Reynolds number Re = 185 have been selected for the physical validation. Here we focus on physical parameters such as hydrodynamic coefficients.

Static cylinder in a cross-flow

We consider in this subsection the classic test case of a static circular cylinder in a cross-flow. The flow is characterized by the Reynolds number Re = UD ν where U is the incoming velocity, D the cylinder diameter and ν the kinematic viscosity of the fluid. The simulations are performed in the unsteady laminar regime with Re = 100. The computational domain Ω, illustrated in the Fig. 17(a), corresponds to a square of length L with a centered immersed cylinder. The boundaries of the computational domain ∂Ω must be located sufficiently far enough to reduce the impact of boundary conditions on vortex development behind the cylinder. Thus, the ratio L/D is set to 60. Symmetry conditions are prescribed on Γ S , imposed velocity at the inlet and imposed pressure at the outlet. The hydrodynamic coefficients considered in this study are the drag C D = F x 0.5ρU 2 D and the lift C L = F y 0.5ρU 2 D with (F x , F y ) t the fluid mechanical constraint integrated on Σ. The results are compared quantitatively in terms of these hydrodynamic coefficients and confront with Introïni & al. [39] that developed a penalized direct forcing method, quite similar to ours, of order 2 in space. Three uniform Cartesian grids are set, ranging from D/∆x = 12.5 to D/∆x = 50 with ∆x the size of the Cartesian cells. Qualitatively, both the standard and the RG formulations of the (base or linear) DF method have the same trend when refining the mesh. That is coherent with the way the RG formulation is implemented. Quantitatively, whatever formulation is considered (standard or RG linear), the time-averaged drag and the drag and lift amplitudes are well calculated with less than 2%-5% of relative deviation for the finest mesh.

Given the results obtained for this test case, it can be emphasized that the RG formulation does not alter the hydrodynamic coefficients values. In fact, the RG linear DF provides even better results than the standard linear DF.

Oscillating cylinder in cross-flow

The effect of the cylinder motion on the hydrodynamic coefficient values is investigated on this test case. We expect that the RG formulation should improve the results in comparison with the original one as it cut off the force oscillations. This test case consists in one of the numerical experiments conducted by Guilmineau and Queutey [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF]: a circular cylinder with an imposed harmonic motion. As for the static cylinder, the flow is characterized by the Reynolds number Re = UD ν . The simulations are performed in the unsteady laminar regime with Re = 185. The boundary conditions and the computational domain are the same as for the static cylinder test case.

The imposed sinusoidal motion is defined by: x c (t) = x c (0) + A e cos(2π f e t) y c (t) = y c (0) where (x c (0), y c (0)) t are the initial coordinates of the cylinder, A e the amplitude of the oscillation and f e its frequency.

Simulations are performed for ratios of A e /D = 0.2 and f e / f 0 = 1.1 where f 0 is the natural shedding frequency from the stationary cylinder at Re = 185. From [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF], we set S t 0 = f 0 D U ≈ 0.19 with U the incoming velocity. Then our reference Strouhal number is S t re f = 1.1 S t 0 ≈ 0.209. Here again, three uniform Cartesian grids are set characterized by: D/∆x = 12.5, D/∆x = 25 and D/∆x = 50. Figs 19(a) and 19(b) give an illustration of the flow physic, displaying the stream lines and the vorticity around the circular cylinder when it reaches its upper position. The results are compared quantitatively in terms of hydrodynamic coefficients and confront with the results from Guilmineau and Queutey [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF]. In this test case, the time-averaged drag coefficient C D and the standard deviations C D,σ and C L,σ of the drag and lift are considered.

Fig. 20 does a comparison between the standard and the RG methods equipped with the base model in term of force oscillations for the coarsest space-step D/∆x = 12.5. It can be clearly noticed that the RG formulation decreases drastically the force oscillations without altering the physics. This is qualitatively confirmed by a comparison between the drag and lift coefficients over a few periods of time obtained by the RG-linear DF method and by Guilmineau and Queutey, cf Fig. 21. Stand Tab. 2 shows the values of the hydrodynamic coefficients computed with the linear model and a grid spacing ratio of D/∆x = 25. The results for all the mesh sizes and both the formulations are displayed on Fig. 22. As a whole, for a given hydrodynamic coefficient, all the computations converge in space toward the same value as expected. The standard and RG DF formulations retain the same behavior when decreasing the grid size. We conclude that the RG formulation is able to reduce the SFO while providing good results in comparison with the literature. The relative deviation of the hydrodynamic coefficients remains lower than 3% (resp. 9%) for C D and C L,σ (resp. C D,σ ). Figure 22: Unsteady flow around a cylinder in a harmonic motion [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF]: (a) time-averaged drag, (b) drag standard deviation and lift standard deviation for both the standard and the RG formulations of the DF method (base and linear models) in comparison with results from Guilmineau and Queutey [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF].

Summary

The main goal of this article was twofold. First, we have driven a theoretical analysis of the dependency to space and time steps of the SFO magnitude observed in the simulations of moving boundary problems with sharp-interface IBM. The SFO decreases by decreasing the grid spacing and increasing the time-step. Second, we have presented a regularized formulation of the DF method, capable of cutting off the SFO. As the standard DF, it consists in solving the Navier-Stokes equations on a Cartesian grid by adding a direct forcing term that mimicks the presence of the solid. But here, the forcing term is regularized through an appropriated weighting. Contrary to the simple solid volume fraction weighting, the expresion of the regularizing function of the forcing term is carefully chosen to smooth the fluid temporal discontinuity through the fluid-solid interface leading to a dramatic reduction of the SFO without altering the physical phenomena at stake. Concerning the design of the immersed boundary condition for the velocity, the RG formulation introduces an explicit-in-time weighting between the values resulting from the interpolation scheme and from the flow-solver numerical scheme. This velocity mixing is explicit: that is no iteration is run to get an implicit mixing between the two velocities as in [START_REF] Luo | A hybrid formulation to suppress the numerical oscillation caused by immersed moving boundaries[END_REF][START_REF] Luo | On numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries[END_REF], leaving the native numerical scheme of a given code unchanged. Here, two interpolation schemes have been considered: the direct solid velocity affectation and the linear interpolation from [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF].

To assess the spatial and time dependencies of the oscillation amplitude, an analysis has been conducted to predict the variation of the forcing term for the standard and RG DF methods with a first-order interpolation model. Concerning the standard DF method, the SFO behavior is found in O( ∆x 2 ∆t ) + O(∆x) due to the DCs that bring the main contribution. For small time-steps, this is in good agreement with the literature: O( ∆x 2 ∆t ) [START_REF] Lee | Sources of spurious force oscillations from an immersed boundary method for moving-body problems[END_REF] or O( ∆x d ∆t ) [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. Concerning the RG DF method, the SFO behavior is found in O(∆t) + O(∆x). As the term O( ∆x 2 ∆t ) is replaced by O(∆t), it comforts the fact that the RG formulation decreases the spurious force oscillations. Then, the analysis results have been confronted to the numerical test case of Seo & Mittal [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF], an oscillating cylinder in a fluid at rest with a Reynolds number of about 78, using both the original and the RG approaches, by quantifying the SFO thanks to the pressure-drag 2-δ discontinuity. The numerical spatial and time dependencies of the oscillations are very close to the theoretical ones predicted by our analysis for the standard and the RG DF methods. From a general point of view, the RG approach decreases the oscillations by at least one to two orders of magnitude over a large range of grid spacing and time-step without increasing the computational cost. The results obtained in this matter are quite similar to [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] that used a cut-cell method that is much more complicated to implement.

Through the numerical simulation of the academical test case of the laminar cylindrical Taylor-Couette flow [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF], it results that the quadratic numerical rate of convergence in L 2 norm of the linear DF method is kept by the RG DF method. At last, two numerical experiments have been conducted to assess the physical validity of the RG formulation: a static cylinder subjected to a cross-flow with a Reynolds number of 100 and a cylinder with an imposed sinusoidal motion subjected to a cross-flow of Reynolds number 185. When the solid is in motion, the RG formulation decreases dramatically the force oscillations especially when the time-step and/or the mesh refinement are quite low (i.e. 10 < D/∆x < 20). Generally speaking, considering the hydrodynamic coefficients, the behaviour of the original method tends to be the same as the RG formulation as expected when refining the mesh. It is coherent with their implementations. These results are quite in good agreement with the literature and prove that the RG linear method is capable of computing complex fluid flows around moving boundaries.

Concerning our field of applications, this work constitutes a first step toward more complicated FSI problems in two-phase flows such as the fluid-elastic instability of moving tubes in a tube bundle subjected to a transverse twophase flow. Indeed, the regularized direct forcing approach adopted in this article is easy to couple with ingredients for this kind of simulations as level-set methods using, for instance, the Desprès-Lagoutière scheme for advection [START_REF] Després | Contact Discontinuity Capturing Schemes for Linear Advection and Compressible Gas Dynamics[END_REF] on the fluid density.

Figure 1 :

 1 Figure 1: Schematic representation of the computational domain Ω and the fluid-solid interface Σ = ∂Ω F ∩ ∂Ω S .

Figure 2 :

 2 Figure 2: Staggered arrangement of the unknowns on a 2-D Cartesian cell (i; j) with the CVs of pressure ((a), blue ) and the CVs of velocity components ((b) and (c), red ).

Figure 3 :

 3 Figure 3: Schematic representation of the base and the linear model on a staggered grid. The red diamonds represent the velocities forced to the solid velocity and the green triangles are the ones forced to the interpolated velocity for the linear model.

Figure 4 :

 4 Figure 4: Seo and Mittal's test case [36]: (a) geometry and computational grid (zoom around the obstacle) and (b) physical dimensionless coefficients time histories: C P , C ν , C T and C DF . Computation done with the DF base model, cf. subsection 2.3.1.

Figure 5 :

 5 Figure 5: C P time histories for the Seo and Mittal's test case: (a) ∆t = ∆t 0 = 0.002T 0 for different grid spacing and (b) D/∆x = 16 for different time-steps.

Figure 6 :

 6 Figure 6: Forcing term maximum for the Seo and Mittal's test case (m/s 2 ). DC and FC numbers on the left and right sides.

Figure 7 :

 7 Figure 7: C 2δ P,σ and C 2δ P,Max for the Seo and Mittal's test case versus: (a) grid spacing (∆t = 1 8 ∆t 0 ) or (b) time-steps (D/∆x = 16).

  in 1D. As the forced-free interface Σ base I of the base model coincides with Σ (d base 0 = 0, cf. Fig. 3), this definition leads to the solid domain fraction ξ base = ∆x+2d(x) 2∆x

Figure 8 :

 8 Figure 8: Representation of χ m (x) and ξ m (x) =

Figure 9 :

 9 Figure 9: Schematic representation of the regularized formulation: a) time step n and b) time step n + 1. The fluid velocities located in the free domain are in green and those in the forced domain Ω m I are in red. The forced-free interface Σ m I is a red line and the velocities where their control volumes are occupied by a fraction of the forced domain are in purple.

Figure 10 :

 10 Figure 10: Schematic representation of a dead cell.

Figure 11 :

 11 Figure 11: C P time histories for the Seo and Mittal's test case and the base model: (a) ∆t = ∆t 0 = 0.002T 0 , for different grid spacing and (b) D/∆x = 16 for different time-steps. DF base: standard DF. RG base: regularized DF.

Figure 12 :

 12 Figure 12: |C 2δ P | time histories for the Seo and Mittal's test case and the base model: (a, b) ∆t = ∆t 0 = 0.002T 0 for different grid spacing and (c, d) D/∆x = 16 for different time-steps. RG (resp. DF) denotes the regularized DF method (resp. the standard DF method).

∆∆Figure 13 :

 13 Figure 13: C 2δ P,σ and C 2δ P,Max for the Seo and Mittal's test case versus grid spacing: (a) ∆t = ∆t 0 = 0.002T 0 , (b) ∆t = 1 4 ∆t 0 and (c) ∆t = 1 8 ∆t 0 .

Figure 14 :

 14 Figure 14: C 2δ P,σ and C 2δ P,Max for the Seo and Mittal's test case versus time-steps: (a) D/∆x = 16 and (b) D/∆x = 32.

Figure 15 :

 15 Figure 15: C p 2δMax for the base and the linear interpolation schemes for both the DF and RG formulations.

Figure 16 :

 16 Figure 16: Cylindrical Taylor-Couette problem : (a) computational domain and (b) relative L 2 (Ω f ) norm of the error vs the ratio of the domain length L over the number of elements. The CFL number is set to 1. Comparison of the base model with the interpolation model for both the standard and the regularized direct forcing.

|ω 1 |r 2 1 ν 2 <

 12 to be set to 1. Then, the Taylor number T a = 0.5Re 2 (r 1 +r 2 )(r 2 -r 1 T a c where T a c ≈ 1.712 is the critical Taylor number[START_REF] Guyon | Physical Hydrodynamics[END_REF]. This assumption implies a strictly 2D flow. Fig. 16(b) presents the L 2 (Ω f ) norm 2 defined by the Eq. (

Fig. 17 (Figure 17 :

 1717 Fig. 17(b) gives an illustration of the type of C D and C L signals obtained by the simulations with D/∆x = 12.5 and the standard base-model DF method. The streamlines and the vorticity are displayed on Figs 17(c) and 17(d).

Figure 18 :

 18 Figure 18: Static cylinder in a cross-flow with Re = 100: (a) time-averaged drag, (b) drag amplitude and (c) lift amplitude for both the standard and the RG formulations. Comparison with the results of Introïni & al. [39] (L/D = 60 and D/∆x = 50).

Figure 19 :

 19 Figure 19: Unsteady flow around a cylinder in a harmonic motion: (a) stream lines and (b) vorticity for the RG linear formulation when the cylinder is at its upper position. The stream lines are colored in function of the velocity along the x axis value. The vorticy red lines are negative and the blue one are positive.

Figure 20 :

 20 Figure 20: Time history of the drag and lift coefficients over a period of time for (a) the standard linear and (b) the RG linear formulations with a computational domain of L/D = 60 and a grid spacing ratio of D/∆x = 12.5; T = 1/ f e .

Figure 21 :

 21 Figure 21: Comparison of the drag and lift coefficients over a few periods of time: (a) picture from the article of Guilmineau and Queutey [47] and (b) results from the RG-linear formulation with a computational domain of L/D = 60 and a grid spacing ratio of D/∆x = 25. T = 1/ f e .

  The maximum and the standard deviation values of C 2δ P , noted C 2δ P,Max and C 2δ P,σ , are calculated over a period of time. Let us notice than C 2δ P,σ is bounded by C 2δ P,Max . Using the same numerical parameters as in Section 3.1, the results are plotted in Figs 7(a) and 7(b) as a function of grid spacing and time-step size. A power-law fitting gives:

	C 2δ P,Max =	∆x 1.9532 ∆t 0.9018 and C 2δ P,σ =	∆x 1.5503 ∆t 0.472 .

Table 1 :

 1 Hydrodynamic coefficients associated with the problem of unsteady flow around a static cylinder, Re = 100 and D/∆x = 25.

			3.1	1.3777	2.3	1.347	1.317 -1.392
	C D	0.01013	12.6	0.009512	5.7	0.009	0.009 -0.012
	C L	0.3429	5.2	0.3338	2.4	0.326	0.303 -0.349
	S t	0.165	< 0.2	0.165	< 0.3	0.165	0.164 -0.172

Table 2 :

 2 Hydrodynamic coefficients associated with the problem of unsteady flow around a cylinder in a harmonic motion[START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF], Re = 185 and D/∆x = 25. Our reference Strouhal number is S t re f = 1.1 S t 0 ≈ 0.209.

		. linear Relative deviation RG linear Relative deviation Guilmineau & al. [47]
			vs [47] (%)		vs [47] (%)
			L/D = 60			L/D = 50
	C D	1.415	0.4	1.366	3.8	1.42
	C D,σ	0.134	10	0.119	20	0.149
	C L,σ	0.844	5.9	0.840	6.4	0.897
	S t	0.203	2.9 vs S t re f	0.206	1.4 vs S t re f	0.214

Acknowledgements References