C Introïni

M Belliard

C Fournier

A Second Order Penalized Direct Forcing for Hybrid Cartesian/Immersed Boundary Flow Simulations

Keywords: immersed boundaries, penalized direct forcing, projection scheme, Cartesian 1 grids, interpolation scheme

In this paper, we propose a second order penalized direct forcing method to deal with fluidstructure interaction problems involving complex static or time-varying geometries. As this work constitutes a first step toward more complicated problems, our developments are restricted to Dirichlet boundary condition in purely hydraulic context. The proposed method belongs to the class of immersed boundary techniques and consists in immersing the physical domain in a Cartesian fictitious one of simpler geometry on fixed grids. A penalized forcing term is added to the momentum equation to take the boundary conditions around/inside the obstacles into account. This approach avoids the tedious task of re-meshing and allows us to use fast and accurate numerical schemes. In contrary, as the immersed boundary is described by a set of Lagrangian points that does not generally coincide with those of the Eulerian grid, numerical procedures are required to reconstruct the velocity field near the immersed boundary. Here, we develop a second order linear interpolation scheme and we compare it to a simpler model of order one. As far as the governing equations are concerned, we use a particular fractional-step method in which the penalized forcing term is distributed both in prediction and correction equations. The accuracy of the proposed method is assessed through 2-D numerical experiments involving static and rotating solids. We show in particular that the numerical rate of convergence of our method is quasi-quadratic.

Introduction

Fluid flow with heat and mass transfer around complex stationary or moving geometries (solid or flexible) appears in a large number of situations of practical interest including biological fluid mechanics (blood flow in human heart for instance) or in life-science context (the fish-like swimming e.g.). Fluid-structure interaction problems are also of importance in many engineering applications, as for example, to design industrial heat exchangers, aerospace vehicles or in nuclear safety context. In this latter case, the vitrification process for the radioactive waste storage is an example. In this process, a viscous multiphase multicomponent flow at high temperature (gas bubbles and molten glass incorporating the ultimate waste) interacts with both static (e.g. the vessel structure, the apparatus of measurement, . . . ) and moving (e.g. the mechanical stirrer) bodies of more or less complex geometries.

The numerical treatment of these kinds of problem appears to be a challenging task because of time-varying geometries, often combined with complex flow regimes. To tackle numerically these complex problems, the well-known body-fitted approach is usually followed.

Such an approach consists in discretizing the governing equations on a non-structured mesh for which the boundaries of the computational domain lie on those of the physical domain.

Thereby, boundary conditions are directly (and so, exactly) imposed on the physical domain boundary. However, the main drawback of the body-fitted like techniques lies in their lack of ability to handle complex industrial problems involving moving bodies which require the development of specific numerical schemes to deal with the difficult issue of re-gridding.

Another approach consists in using non-boundary conforming techniques in which the physical domain is immersed in a fixed fictitious one of simpler geometry on a Cartesian grid. Such techniques allow us to use efficient, fast and accurate numerical methods avoiding the tedious task of the re-meshing caused by time-varying geometries. In contrast, as the immersed boundaries are described by a set of Lagrangian points (or the zero of a level-set function) that do not generally coincide with those of the Eulerian grid, numerical methods have to account the immersed boundary conditions at their right places. The non-boundary conforming techniques proposed in the literature may be classified into two categories.

The first category, including for instance Cartesian methods (e.g. [START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF][START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF]), the Immersed Interface Method (IIM, [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF]) or the Jump Embedded Boundary Condition method (JEBC, [START_REF] Ramière | A general fictitious domain method with immersed jumps and non-conforming structured mesh[END_REF]), mimics the presence of embedded geometries by modifying the numerical scheme in the immediate vicinity of the immersed boundary or interface. The two latter methods introduce jump conditions across the interface in the solve of the partial differential equations. Such an approach leads to a sharp representation of the immersed interface but, for the Cartesian method, extending it to three-dimensional problems may appear to be a challenging task, particularly regarding the coding logistic.

In the second category (rather than locally modifying the numerical scheme) a supplementary term, referred to as the forcing term, is added to the governing equations. This class of non-boundary conforming techniques dates back to Peskin's works in which an Immersed Boundary Method (IBM) has been developed to numerically simulate blood flows in a human heart [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF]. In this case, the immersed boundaries correspond to muscular heart walls and extra forces acting in these boundaries are modeled by a vectorial forcing term added to the continuous Navier-Stokes equations. A Lagrangian coordinate system is employed to track the interface and to calculate the vectorial forcing term. The IB method has been successfully applied to problems with elastic geometries but, in the rigid limit, it generally leads to very stiff problems. Moreover, in order to ensure the stability of the numerical scheme, the forcing term based on a Dirac delta function must be smeared over a stencil of few Cartesian nodes. Following the ideas introduced by Peskin, several IB-like methods with different forcing terms (or forcing strategies) have been proposed in the literature. In [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF], Goldstein & al. propose the Feedback Forcing (FF) method in which the forcing term can be viewed as a force density that brings the fluid velocity to zero near the immersed boundary.

Similarly to what is done in [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF], the numerical scheme used in [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF] requires a spreading of the forcing term over the interface. Moreover, the FF method suffers from the fact that the forcing term highly depends on flow properties. Whether the Peskin's IB method or the FF method, their application to flows at high Reynolds number is limited by the spreading of the forcing term over the immersed boundary. In this case, local mesh refinement techniques can be a solution [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF]. An alternative approach to the aforementioned techniques, referred to as Direct Forcing (DF) method, has been proposed by Mohd-Yusof [START_REF] Mohd-Yusof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF] and then adapted by Fadlun & al. [9]. This immersed boundary technique consists in directly applying the desired boundary conditions on Cartesian nodes close to the interface leading to a quasi sharp representation of the interface (through one cell layer). In that sense, using the terminology employed by Gilmanov & al. in [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF], the DF method may be referred to as a Hybrid Cartesian/Immersed Boundary (HCIB) approach and may be conceptually related to the IIM. Moreover, one of the interest of the DF method is that the forcing term can be easily computed and it does not depend on the flow properties. Therefore, the stability of the numerical scheme is not affected. However, the accuracy of the DF method is partially dependent on the numerical scheme because the calculation of the forcing term is in particular based on the discretized form of the governing equations. Since its development by Mohd-Yusof [START_REF] Mohd-Yusof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF], the DF method has gained in popularity and has been successfully applied to various fluid-structure interaction problems (e.g. [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF][START_REF] Peller | High-order stable interpolations for immersed boundary methods[END_REF][START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[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][START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF][START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF]) or turbulent flow simulations (e.g. [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF]17]) using mesh refinement or mesh stretching techniques. It is also worth to mention the immersed boundary method of Pinelli & al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume solvers[END_REF], that has roots in both IBM and DF methods, and which is suitable for general grid systems including curvilinear ones.

More recently, Belliard & Fournier [19] have proposed a variant of HCIB techniques, called Penalized Direct Forcing (PDF) method, that combines both the basic features of the DF method and those of L 2 -penalty methods (e.g. [20]). Links can be found with the works of Sarthou & al. [21] and those of Bergmann & Iollo [22]. As for the DF method, the unknowns are locally enforced on the grid nodes nearest the immersed interface. However, the PDF method appears to be a more versatile approach than the DF method because the forcing term expresses as a L 2 -penalty term that is independent on the discrete governing equations.

In the present paper, after introducing the discretization of the Navier-Stokes governing equations in the Section 2, the PDF algorithm is detailed in the Section 3, including a specific treatment of the pressure near the immersed boundaries. The PDF method itself is presented in Section 3.1. Interesting for practical purposes, the non-boundary conforming approaches are often coupled with fractional-step schemes (e.g. [23,24]) but, as emphasized by Ikeno & Kajishima [25] or Taira & Colonius [START_REF] Taira | The immersed boundary method: a projection approach[END_REF], most of them take account of the forcing term only in the prediction equation leading to inconsistent schemes. Here, an original fractional-step scheme leading to a consistent (in the sense of [25]) PDF method is developed and presented in Section 3.2. Homogeneous Neumann IBCs for the pressure are recovered through a particular treatment of the pressure equation coefficients near the immersed interface.

Whatever the non-boundary conforming method involved, an important issue concerns the reconstruction of the velocity field close to the immersed boundary and the accuracy of the numerical method developed for this purpose. These points have received a particular attention in the literature with, most of the time, the development of interpolation schemes.

Most of the classical approaches consist in interpolating or extrapolating the velocity field in a preferred direction (e.g. [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF][START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF]). In this work, we have developed an original robust interpolation scheme, second-order accurate in space, that is not guided by particular direction. It relies mainly on an averaged reconstruction of the velocity gradient near the IB and on an approximate projection operator onto the IB. Without loss of generality, we restrict our presentation to Dirichlet's IBCs for the velocity1 . This is the object of Section 4.

Finally, in Section 5, some 2-D numerical experiments are performed for steady and unsteady incompressible laminar flows at very moderate Reynolds number (up to 100) around/ between static and rotating solids to assess the validity, accuracy and the ability of the proposed method for both uniform and analytically velocity prescribed at IBs. We show in particular that the numerical rate of convergence is (quasi-)quadratic for all studied cases.

To highlight the ability of our method to deal with 3-D moving geometries, we also present an illustration of flows induced by a stirrer .

Governing equations and numerical method

The section is devoted to the numerical method. It is structured in two parts. The first part focuses on the governing equations whereas the numerical scheme is the object of the second part.

Governing equations

The governing equations used to describe unsteady incompressible flows are given by:

∂u ∂t + ∇ • (u ⊗ u) + ∇P -ν∇ 2 u = f in Ω (1a) ∇ • u = 0 in Ω (1b) u = u D on ∂Ω and u(t 0 ) given in Ω (1c)
where Ω denotes the computational domain, ∂Ω its boundary, u the solenoidal velocity and ν the kinematic viscosity. Here-above, P is the total pressure defined by:

ρ∇P = ∇p -ρg (2)
where p is the hydrodynamic pressure, ρ the constant density and g the gravity force. For the sake of clarity, we assume full Dirichlet boundary conditions and we consider the volume force f = 0 in the sequel of the paper.

Numerical scheme

Here, we focus on the numerical scheme used in this work. On the first hand, we present the space and time discretizations. On the second hand, we discuss about the fractional-step method used to solve the Navier-Stokes equations.

Time and space discretizations

The time marching of the velocity is performed by means of a variant of the degenerate fourth-order explicit Runge-Kutta scheme of Williamson [START_REF] Williamson | Low-storage Runge-Kutta schemes[END_REF]. Given an initial condition u 0 = u(t 0 ) and ∂u ∂t = f (u), the new velocity u n+1 = u(t 0 + (n + 1)∆t) is obtained by:

∀ n ∈ N, u n+1 = u n + β 1 q 1 + β 2 q 2 + β 3 q 3 (3)
where

β 1 = 1 2 and q 1 = ∆tf (u n ) , β 2 = 1 and q 2 = ∆tf u n + 1 2 q 1 - 1 2 q 1 , β 3 = 1 6 and q 3 = ∆tf u n + 1 2 q 1 + q 2 -2q 2 .
This scheme is degenerated because it only uses two storage vectors and three computational steps like the third-order Runge-Kutta scheme. In our variant, at each step k of the Runge-Kutta scheme, we get f (u n +. . .) by solving the couple (u, P ) using a fractional-step method (cf. Section 2.2.2) for which the temporal discretization is based on a semi-implicit scheme which relies on an explicit discretization of the convection term and an implicit discretization of the diffusion.

As far as the space discretization is concerned, it is based on a finite volume approximation with a staggered grid arrangement of the primitive 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.

[Figure 1 about here.] Fig. 1 presents an example of the staggered arrangement of the unknowns on a 2-D Cartesian cell Ω i,j . In this figure, Ω i+1/2,j and Ω i,j+1/2 are respectively the control volumes of the components u i,j and v i,j of the velocity field. Note that the control volume of the pressure P i,j coincide with the cell Ω i,j . In this frame, the governing equations Eqs. (1) are integrated over each control volume ensuring the conservation of mass and momentum balances. Moreover, the convection and diffusion terms are respectively approached by the QUICK and the centered schemes [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. Now, denoting by ⋆ h a discrete space operator associated to ⋆, at each step k = 1, . . . , 3 of the Runge-Kutta scheme (3), the discrete form of the governing equation Eqs. (1) reads:

f (u k-1 ) := u ⋆ -u k-1 ∆t = -∇ h • u k-1 ⊗ u k-1 -∇ h P ⋆ + ν∇ 2 h u ⋆ in Ω (4a) ∇ h • u ⋆ = 0 in Ω (4b) u ⋆ = u D on ∂Ω (4c)
where

u k = u n + k i=1
β i q i , u 0 = u n and (u n+1 , P n+1 ) = (u 3 , P ⋆ ). Here, (u ⋆ , P ⋆ ) stands for a solenoidal velocity and a pressure obtained by a projection method as described in the following part.

A fractional step algorithm to solve the Navier-Stokes equations

The fractional-step method or projection method has been introduced by Chorin and

Temam in 1968 for incompressible flows [23,24]. On the basis of their work, many variants have been proposed, as for example, the incremental projection method for incompressible flows [START_REF] Shen | On Error Estimates of some Higher Order Projection and Penalty-Projection Methods for Navier-Stokes Equations[END_REF], the projection schemes for dilatable or barotropic fluids [START_REF] Jobelin | Une méthode de pénalité-projection pour les écoulements dilatables[END_REF][START_REF] Th | An unconditionally stable pressure correction scheme for the compressible barotropic Navier-Stokes equations[END_REF] or, more recently, the novel fractional time stepping technique massively parallel for incompressible Navier-Stokes equations developed by Guermond & Minev [START_REF] Guermond | A new class of massively parallel direction splitting for the incompressible Navier-Stokes equations[END_REF], to cite among others. We refer the reader to [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for a recent review of these methods.

In this paper, we use a non-incremental fractional-step scheme. The first step consists in solving a predicted velocity u without pressure gradient term as follows2 :

u -u n ∆t + ∇ h • (u n ⊗ u n ) -ν∇ 2 h u = 0 in Ω (5a) u = u D on ∂Ω. (5b)
Then, the second step corresponds to a correction stage which consists in computing a new pressure P n+1 and recovering a new solenoidal velocity u n+1 . By assuming that ∇ 2 h u ∼ ∇ 2 h u n+1 , this step reads:

u n+1 -u ∆t = -∇ h P n+1 in Ω (6a) ∇ h • u n+1 = 0 in Ω (6b) n • u n+1 = n • u = n • u D on ∂Ω. (6c)
Using Eq. (6b) and Eq. (6c), one can write:

∇ h • ∆t∇ h P n+1 = ∇ h • u in Ω (7a) n • ∆t∇ h P n+1 = 0 on ∂Ω. (7b)
Finally, the correction Eqs. ( 6) allows us to compute the new velocity u n+1 as follows:

u n+1 = u -∆t∇ h P n+1 in Ω. ( 8
)

Immersed boundary method

The root concept of immersed boundary approaches is depicted in Fig. approach, we refer the interested reader to [START_REF] Unverdi | A Front-Tracking method for viscous, incompressible, multi-fluid flows[END_REF][START_REF] Tryggvason | A Front-Tracking method for the computations of multiphase flow[END_REF], to cite among others.

[Figure 2 about here.]

The remain of this section is structured as follows. In the first part, we introduce the concept of the penalized direct forcing and, in the second part, the consistent fractionalstep algorithm used with the PDF method. The interpolation scheme developed in order to reconstruct accurately the velocity field near the immersed boundary is presented in Section 4. Let us remark that, even if the algorithm is illustrated on fixed rigid bodies on uniform Cartesian grids, it also works on moving bodies and on non-uniform grids.

Furthermore, it could be extended to computational grids that can be mapped into Cartesian ones or to deformable bodies as in [22].

Penalized direct forcing

In the frame of immersed boundary-like methods [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF][START_REF] Mohd-Yusof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF], the presence of embedded time-varying geometries Ω s in the computational domain Ω is taken into account by a supplementary source term F added to the Navier-Stokes equations given by Eqs. [START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF]. This leads to the following equations:

∂u ∂t + ∇ • (u ⊗ u) + ∇P -ν∇ 2 u = F in Ω (9a) ∇ • u = 0 in Ω. ( 9b
)
The term F, here after referred to as penalized direct forcing term, is defined by [19]:

F = αχ s η (u imp -u) with α > 0, 0 < η ≪ 1 and χ s → [0, 1] (10)
where χ s is the characteristic function or the volume ratio function of the obstacle [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF], η a penalty coefficient, 0 < α a real parameter and u imp the imposed fluid velocity around/inside the obstacle Ω s .

The imposed fluid velocity u imp depends not only on the solid velocity u s that may be uniform (i.e. constant in time and space) or analytic (i.e. u s := u s (x, t)) but also on the fluid velocity u depending on the requested accuracy. This dependency leads to a change of the mass matrix when solving the equation Eq. (9a). But in practice, such a modification appears to be a complex and tedious task as suggested by the relationship Eq. ( 26) of Section 4. In this work, rather than tackle directly this full difficulty appearing in equation Eq. (9a), our objective is to calculate the velocity u n+1 using an explicit approximation of u imp in the forcing term. This is achieved by solving first the following free-obstacle explicit N.-S.-like equations3 :

ũ⋆ -u n ∆t ⋆ + ∇ h • (u n ⊗ u n ) -ν∇ 2 h u n = 0 in Ω. (11)
Then, the resulting provisional velocity ũ⋆ is used to calculate an approximate imposed fluid velocity u n+1 imp that depends on u n+1 s and ũ⋆ . A similar strategy can be found in [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume solvers[END_REF].

Remark 1. Eq. ( 10) can be viewed as an implicit limit version (∆t = η/α ≪ 1) of the direct forcing expression

F DF = χ s ( u imp -u k-1 ∆t -f (u k-1
)) [START_REF] Mohd-Yusof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF] applied to the Eq. (4a):

F DF = χ s ∆t {u imp -(u k-1 + ∆tf (u k-1 ))}. ( 12
)
Remark 2. As far as the penalty coefficient η is concerned, Angot & al. [20] have shown that the solution of the PDE system given by Eqs. The penalized forcing term F is only applied on Cartesian nodes near or inside the immersed boundary (i.e. χ s > 0; at the discrete level χ s is taken as an extension of the cell volume ratio function of the obstacle strictly positive for the penalized faces). On these nodes, thanks to the penalty coefficient η -1 , the momentum equation Eq. (9a) boils down to the following relationship:

u(x) = u imp (x), for x ∈ Ω such that χ s (x) > 0. ( 13
)
Note also that the classical incompressible Navier-Stokes equations are recovered in the (computed) fluid domain Ω f when χ s is set to zero in Eq. [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF].

At this point of the analysis, it remains to define the imposed fluid velocity u imp at Cartesian nodes nearest Σ h . Regarding those that are located into the solid domain Ω s (depicted by ⋄ in Fig. 3), the imposed velocity u imp is exactly set to be the obstacle velocity u s . This can be express as:

u imp (x) = u s (x), for x ∈ Ω such that χ s (x) > 0 and d(x) < 0 ( 14
)
where d(x) denotes a signed distance to the immersed boundary Σ h4 .

[Figure 3 about here.]

As far as the other nodes are concerned, namely the fluid nodes belonging to the interfacial region with a positive signed distance (i.e. χ s (x) > 0 and d(x) > 0), a numerical procedure is required to reconstruct accurately the velocity field. This is the object of Section 4.

A consistent fractional-step method

In the literature, most of the fractional-step schemes, employed in the frame of immersed boundary methods, do not take account of the forcing term in the correction step. As emphasized by Ikeno and Kajishima [25], such an approach leads to inconsistent schemes in the sense that the immersed boundary conditions are well satisfied by the predicted velocities in the prediction stage but not by the new velocities at the end of the projection stage. See also [START_REF] Taira | The immersed boundary method: a projection approach[END_REF] for a projection approach of the immersed boundary methods.

In this section, we propose a modification of the fractional-step method presented in Section 2 in order to obtain a consistent scheme in the sense defined in [25].

The new feature of our algorithm [19] is that the forcing term is distributed both in the prediction and the correction stages of the projection. Under these circumstances, the first step of our scheme consists in solving a predicted velocity u without the pressure gradient term as follows (setting α = 1/∆t):

u -u n ∆t + ∇ h • (u n ⊗ u n ) -ν∇ 2 h u = χ s η∆t (u n+1 imp -u) in Ω. (15)
Here, as in the standard fractional-step algorithm, u respects the following boundary conditions on the obstacles:

u = u D on Σ h . (16)
The second step of our modified algorithm corresponds to the correction stage to recover the pressure P n+1 and the solenoidal velocity u n+1 , but including the remaining part of the forcing term:

u n+1 -u ∆t = -∇ h P n+1 + χ s η∆t ( u -u n+1 ) in Ω (17a) ∇ h • u n+1 = χ s ∇ h • u n+1 imp in Ω. ( 17b
)
The equation Eq. (17a) can be simplified as:

ρu n+1 -u ∆t = -∇ h P n+1 in Ω. ( 18
)
Here, it is interesting to remark that the equation Eq. ( 18) is similar to the classical correction equation Eq. (6a) with ρ = (1 + χ s /η). This justifies our choice α = 1/∆t. Now, using Eq. (17b), we obtain5 :

∇ h • ∆t ρ ∇ h P n+1 = ∇ h • u -χ s ∇ h • u n+1 imp in Ω. ( 19
)
In the liquid Ω f (i.e. χ s = 0), we have ρ = 1. But in the solid (i.e. χ s > 0), it becomes ρ(η) ∼ O η -1 . Therefore, 1/ρ(η) can be viewed as an effective coefficient which satisfies:

1 ρ(η) ∼ O (η) ≪ 1. ( 20
)
By doing this, we enforce the local homogeneous Neumann boundary conditions of Eq. ( 19)

on obstacles in a natural way. The use of small values for the diffusion coefficients to simulate immersed homogeneous Neumann BCs can also be found in [START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF].

In Eq. ( 19), in order to avoid numerical difficulties caused by these low diffusion coefficients inside the immersed boundary (when all the velocity degrees of freedom are penalized), we may add a pressure L 2 -penalty term leading to the following modified pressure correction equation:

∇ h • ∆t ρ ∇ h P n+1 + χ s η P 0 -P n+1 = ∇ h • u -χ s ∇ h • u n+1 imp in Ω. ( 21
)
where P 0 is a prescribed pressure correction. This only affects the cells totally included in the obstacle for which the velocities and the pressures are imposed but not the consistency of the scheme in the fluid sub-domain. In practice, P 0 is set to zero (cf. Sec. 5).

Remark 3. To clarify this point, let us consider as example an 1D case with the cell i in fluid and the cell i + 1 (cut by Σ) in the solid. The face i + 1/2 coincides with Σ h and, following Remark 2, we have:

u Σ h = u n+1 imp Σ h + O(η ζ ) with ζ = 3/4 in the approximated-solid
region. On one hand, Eq. ( 19) and Eq. ( 21) are similar for the fluid cell i and give us:

u i-1/2 -∆t ρ (∇ h P n+1 ) i-1/2 = u n+1 imp Σ h + O(η ζ ) -∆t ρ O(η)(∇ h P n+1 ) Σ h . This does not break the imposed boundary condition on Σ h if (∇ h P n+1 ) Σ h is O(η ζ-1
) or higher order.

On other hand, Eq. ( 21) applied to the solid cell i + 1 gives us:

∆t ρ O(η)(∇ h P n+1 ) i+3/2 - ∆t ρ O(η)(∇ h P n+1 ) Σ h + 1 η P 0 -P n+1 i+1 = O(η ζ ) leading to P n+1 i+1 = P 0 +O(η 1+ζ
) and (∇ h P n+1 ) Σ h as well as (∇ h P n+1 ) i+3/2 = O(η ζ-1 ) or higher orders.

Finally, the correction equation ( 18) allows us to compute the new velocity u n+1 as follows:

u n+1 = u - ∆t ρ ∇ h P n+1 in Ω. ( 22
)
As a whole, we get a consistent fractional-step scheme in the sense that the immersed Dirichlet boundary condition is well satisfied on the obstacles by the new velocity u n+1 , see Eqs. ( 16) and ( 22):

n • u n+1 = n • u = n • u n+1 imp on Σ h ( 23
)
as well as the immersed Neumann boundary condition by the pressure, see Eq. ( 20):

∆t ρ ∇ h P n+1 • n = O (η) on Σ h . (24)

Reconstruction of the velocity field : interpolation schemes

The immersed boundaries are described by a set of Lagrangian points which do not generally coincide with the nodes of the Eulerian mesh. Under these circumstances, the velocity field must be reconstructed near the immersed boundary in order to take account of the immersed boundary conditions. This issue, which constitutes one of the major difficulty of the immersed boundary methods, has received a particular attention in the literature.

[Figure 4 about here.]

The simplest approach consists in prescribing directly the velocity of the obstacle u s on the Cartesian nodes closest to the immersed boundary Σ h without any interpolation scheme [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF].

This numerical procedure, referred as the base model, is first-order accurate in space and reads:

u imp (x) := u s (x) + O (h) , for x ∈ Ω such that χ s (x) > 0 and d(x) > 0 ( 25
)
where h is the size of a Cartesian cell. As illustrated in Fig. 4, for instance with a staggered grid, the base model leads to a step-wise description of the immersed interface Σ h and thus, to a larger obstacle.

Improved numerical algorithms have been proposed in the literature to reconstruct the velocity field in the interfacial region. Among them, the most widely employed to improve the accuracy of the solution close to Σ h are those based on an interpolation or extrapolation procedure (e.g. [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF][START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF][START_REF] Peller | High-order stable interpolations for immersed boundary methods[END_REF][START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[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]25]). Roughly speaking, the interpolation or extrapolation procedures involve solid and fluid velocity contributions to calculate the velocities at the forcing nodes.

The solid contribution corresponds to a node located on the immersed boundary Σ h and is required when the solid velocity is analytic (as u s depends on the position considered on Σ h ). As mentioned in the introduction of this paper, a classical approach consists in interpolating or extrapolating along a specific direction. For instance, in [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF], the solid contribution corresponds to the orthogonal projection of the forcing node onto Σ h . And In this paper, we develop a linear interpolation scheme, second-order accurate in space, in which the fluid contribution is built following similar ideas used in [25] through an averaged reconstruction of the velocity gradient near Σ h . The solid contribution is determined by means of a minimization problem. The scheme consists in an estimation of the imposed fluid velocity u imp at the forcing or penalized node x according to:

u imp (x) := u s (Π Σ (x)) + d(x) N N p=1 u(x p ) -u s (Π Σ (x p )) d(x p ) + O h 2 for x ∈ Ω such that χ s (x) > 0 and d(x) > 0 ( 26
)
where x p = x ± h is the p-th fluid node in the immediate vicinity of the penalized node

x. Fig. 5 shows a schematic representation of the penalized nodes x (red crosses ×) and their associated fluid nodes x p (blue circles •) involved in our interpolation procedure. The relationship of Eq. ( 26) is similar to the one proposed in [25]. It can be viewed as an averaged linear interpolation over all the fluid nodes • located in the immediate vicinity of the penalized node x. Therefore, contrary to the classical approach found in [START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[END_REF] and [START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF],

the local influence of the fluid flow at x is fully taken into account.

[Figure 5 about here.]

In Eq. ( 26), Π Σ (x) denotes the projection of x into the immediate neighborhood V Σ h of the immersed boundary Σ h and corresponds also to the solid contribution of our interpolation scheme. It is defined throughout an algorithm based on the following minimization problem.

Find

z (= Π Σ (x)) ∈ V Σ h such that J(z) = inf y∈V Σ h J(y) (27)
where J(•) is defined by:

∀ y ∈ V Σ h , J(y) = y -x 2 . ( 28
)
To tackle this minimization problem, regardless of the Cartesian node involved x (generic name for forcing nodes or fluid nodes), the first step consists in the partial reconstruction of the immersed boundary in the immediate vicinity of x. To do so, we collect all the Lagrangian facets that belong to the cells intersected by Σ h and located in the immediate vicinity of x. If we consider an interfacial node x (i.e. χ(x) > 0 and d(x) > 0, cf. Fig. 6),

we begin to detect the first level of neighboring cells intersected by Σ h . Then, we identify the "neighbors of neighbors" themselves cut by Σ h and, for each of them, we collect the associated Lagrangian facets. If we consider a fluid node x (i.e. χ(x) = 0), the collecting procedure differs slightly because, in such a case, the neighboring cells of x are always purely fluid cells (cf. Fig. 7). Therefore, the collecting procedure starts with the identification of the interfacial nodes located in the neighboring cells of x. Then, for each interfacial nodes, we identify the neighboring cells cut by Σ h and we collect the associated Lagrangian facets.

[Figure 6 about here.]

[Figure 7 about here.]

After collecting a set of Lagrangian facets, the next step is to determine the equation of the plan passing through each of them. For example, if we consider the i-th Lagrangian facet associated with x, such an equation reads:

d j=1 c ij y j = f i ( 29
)
where d is the spatial dimension of the problem, c ij the coordinates of the unit normal vector at the centroid of the i-th facet and y j the coordinates of one vertex belonging to the i-th facet. Now, by applying Eq. ( 29) on all the facets associated with the Cartesian node x, the minimization problem given by Eq. ( 27) can be expressed as:

Find z (= Π Σ (x)) ∈ V Σ h such that J(z) = inf C•y≈f J(y). (30)
To solve this new minimization problem, we use an Uzawa algorithm which can be summarized as follows.

1. Initialization (k = 0) : we assume λ 0 = 0 2. k-th iteration : by assuming λ k known, we are able to compute ⋆ y k+1 by solving

y k+1 = x - 1 2 C • λ k ⋆ λ k+1 by solving λ k+1 = max λ k + ̺ C • y k+1 -f , 0 with ̺ = ( C 1 C ∞ ) -1
where λ k is the Lagrange multiplier vector at the k-th iteration associated with the constraint

C • y k+1 = f.
Regarding the coding logistic, the computational cost of the Uzawa algorithm depends on the size of the coordinates matrix C. Thus, it depends on the number of Lagrangian facets collected for a Cartesian node x. The collecting procedure, based on a two-level detection of the Cartesian cells crossed by Σ h , contributes to reduce the number of facets without limiting the accuracy of the proposed immersed boundary method (cf. Sec. 5). But, there may be cases in which several Cartesian cells contain a significant number of Lagrangian facets, for instance when the mesh is generated by CAD. To overcome this difficulty, we have implemented a numerical procedure that selects the most representative facets and also detects those that are collinear in order to remove the duplicate. In this study, the most representative facet means the most large. It is worthy to mention that different criteria can be chosen. Our collecting approach also limits the approximation errors in the Uzawa algorithm. Indeed, as the projection of the Cartesian node x is performed onto a set of plans and not onto a particular facet, increasing the number of level (and thus the number of facets) would necessarily increase the number of equations that must be satisfied by the projected node. The proposed algorithm is very robust and z (= Π Σ (x)) ∈ V Σ h converges toward z ∈ Σ h when the space step decreases. Its implementation is a bit less straightforward than in the case of an usual local interpolation algorithm but it allows a great robustness for the complex immersed boundary geometries that can be found in industrial processes.

Numerical experiments

This section is devoted to the numerical validation of the proposed Penalized Direct Forcing method. On the one hand, a Poiseuille flow in an inclined channel, a cylindrical

Couette flow and a steady flow around static/rotating cylinders are considered to assess the validity and the ability of our method both for uniform and analytic Dirichlet's IBCs. All these academic numerical tests have been performed with both the base model (cf. Eq. ( 25))

and the linear interpolation scheme (cf. Eq. ( 26)). Moreover, a grid convergence study has been done in order to obtain the numerical rate of convergence of the method. To do so, we have calculated the relative L 2 (Ω f ) norm ε 2 and the absolute L ∞ (Ω f ) norm ε ∞ of the error of the velocity component u in the fluid domain Ω f as follows:

ε 2 = N f i=0 (u i -u i,ref ) 2 N f i=0 u i,ref 2
and ε ∞ = max 0≤i≤N f (|u i -u i,ref |) (31)
where the lowerscript i denotes the i-th face of the Cartesian grid, N f the total number of faces in the fluid region Ω f , u i the u-velocity component at the center of the i-th face and u i,ref a reference velocity at the same location. This latter term is defined either by the analytic solution in the case of the Poiseuille and Couette flows or by the thirdorder extrapolation of the velocity computed on the finest grid in the case of steady flows around static/rotating cylinders. On the other hand, interestingly for practical purposes, our method coupled with the linear interpolation scheme is also applied to 2-D calculations involving a unsteady laminar flow and to 3-D calculations of the flow induced by a stirrer illustrating the ability to deal with moving geometries.

All the calculations have been performed using the CFD code Trio U, which is a modular (i.e. finite differences or finite volumes) software package for thermal-hydraulic computations, developed at "Commissariat à l' Énergie atomique et aux énergies alternatives" [START_REF] Calvin | A development framework for parallel CFD applications: Trio U project[END_REF].

As a whole, the linear interpolation scheme is about 10%-20% more expensive in CPU time than the base model. Unless specifically mentioned, we solve the arising linear systems from Eqs. ( 15) and ( 19) with the conjugate gradient (CG) method. Considering the DF method, it is enough to precondition the pressure's iterative solver by a SSOR. But for the PDF method, due to badder matrix condition numbers, we need more powerful preconditioners.

We use the diagonal inverse one (low cost) and the algebraic multigrid one for respectively Eqs. ( 15) and (19). This can induce a CPU time over-cost of about 10% respectively to the DF method, but the algebraic multigrid preconditioner is also very useful for the DF method, saving many CG iterations in case of a large number of degrees of freedom.

Poiseuille flow in an inclined channel

Here, the problem under consideration is the well-known test case of a Poiseuille flow in an inclined channel (θ = π/4). We refer to Fig. 8 in which we have shown the geometrical features of the computational domain Ω = Ω s ∪Ω f with the immersed boundary Σ = Ω s ∩Ω f .

[Figure 8 about here.]

Regarding the boundary conditions imposed on Γ s and Γ out , we prescribe a null pressure.

While, we prescribed on Γ in the following parabolic velocity:

u in = (u in , v in ) T = U ∞ -Y 2 , 0 T with Y = y -x √ 2 (32)
where U ∞ = 0.605 m.s -1 is the maximum velocity. Here-above, (X, Y ) is the system of coordinates in the frame associated with the inclined channel. For this numerical test, a no-slip boundary condition is imposed on the immersed boundary Σ (i.e. u s = 0).

[Figure 9 about here.]

Results are the following. Fig. 9 shows the stationary velocity and pressure fields in the computational domain obtained by the proposed penalized direct forcing method with the linear interpolation model. As mentioned previously, we have performed a grid convergence study. Four computational grid sizes have been considered ranging over 25.10 -3 m ≤ h ≤ 2.10 -1 m. Fig. 10 presents the evolution of the L 2 (Ω f ) norm ε 2 and the L ∞ (Ω f ) norm ε ∞ of the error given by Eq. [START_REF] Th | An unconditionally stable pressure correction scheme for the compressible barotropic Navier-Stokes equations[END_REF]. As expected, the graphs of Fig. 10 show that the proposed linear interpolation model leads to a quadratic numerical rate of convergence while the numerical order of the method using the base model is about one.

[Figure 10 about here.]

Here, it is interesting to remark that, in the case of the linear interpolation model, the last three errors in L 2 (Ω f ) and L ∞ (Ω f ) norms are at least one order of magnitude lower than all the ones obtained with the base model. These results are confirmed by the graphs presented in Fig. 11. Indeed, we remark that the velocity profile obtained in the case of the coarse grid with the linear interpolation model is very close to the one computed with the base model on the fine grid. In contrary, the result obtained with the base model on the coarsest grid is very bad. The fact is that the immersed interfaces do not coincide with the grid lines and the base model consists in directly imposing the solid velocity u s on the edges of the cells cut by the immersed interfaces (cf. Eq. ( 25)). This results in an erroneous step-wise description of the walls with a lower flow section. This is no longer the case using the linear interpolation model.

In conclusion, for this test, the results clearly show that the proposed method coupled with the linear interpolation model allows us to obtain an accurate solution at a low computational cost.

[Figure 11 about here.]

Taylor-Couette problem

Here, we focus on a steady flow between two rotating concentric cylinders, also called Taylor-Couette flow. This problem has already been considered in the frame of immersed boundary-like methods (e.g. [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 us to numerically estimate the order of accuracy of the proposed method on problems involving rotating geometries.

Fig. 12 presents the geometrical features of the computational domain Ω = Ω s ∪ Ω f with the immersed boundaries Σ 1 and Σ 2 mimicking the inner-and outer cylinders, respectively. Ω = [0, L] × [0, L] corresponds here to a square where L will be defined later. The inner cylinder rotates clockwise (ω 1 > 0) while the outer cylinder rotates counterclockwise (ω 2 < 0). In overall calculations, we assume that the Reynolds number, defined by Re = |ω 1 |r 2 1 /ν, is set to 1. Such an assumption allows us to write:

Ta < Ta c ( 33
)
where Ta = 3/2 is the Taylor number defined by Ta = 0.5Re 2 (r 1 + r 2 )(r 2 -r 1 ) 3 /r 4 1 and Ta c = 1.712 is the theoretical critical one [START_REF] Guyon | Physical Hydrodynamics[END_REF]. Eq. ( 33) implies that the following simulations are still strictly planar (2-D).

[Figure 12 about here.]

The numerical parameters of the computations are summarized in Tab. 1.

[Table 1

about here.]

Regarding the boundary conditions, we prescribed symmetry conditions on the boundaries of the computational domain Γ and the following analytic velocities on the immersed boundaries Σ 1 and Σ 2 :

∀ i = 1, 2 u s = ω i r i e r on Σ i ( 34
)
where e r is the radial unit vector. The flow is initially at rest and converges toward a steady state. Fig. 13 presents an example of uniform Cartesian grid used to perform the calculations (left) and the steady velocity field (right).

[Figure 13 about here.]

Here again, a grid convergence study has been done in order to estimate the numerical order of the proposed method on problems involving moving bodies. For this problem, five meshes have been considered over the range 25.10 -3 ≤ h/r 1 ≤ 10 -1 where h denotes the step size.

In order to reduce the computational time of these simulations, we have defined the length L of the computational domain Ω as follows:

L = 2 (r 1 + r 2 + Nh) ( 35
)
where N is the number of cells in the outer cylinder (here, N = 3). This strategy is justified because, in the solid, the imposed fluid velocity is equal to the solid velocity (cf. Eq. ( 14)).

As previously done, the linear interpolation model is compared to the base model. Fig. 14 presents the L 2 (Ω f ) norm ε 2 and the L ∞ (Ω f ) norm ε ∞ of the error defined by Eq. [START_REF] Th | An unconditionally stable pressure correction scheme for the compressible barotropic Navier-Stokes equations[END_REF].

As in the case of an uniform imposed velocity (cf. part 5.1), a quasi-linear numerical rate of convergence is obtained with the base model and a quasi-quadratic rate with the linear interpolation scheme. These good results confirm the enhancement of the accuracy as expected with the linear interpolation model. However, it is more difficult to explain why, in this case, the numerical rate of convergence is slightly lower than two (cf. Fig. 14). Such a result might be a consequence of two difficulties associated with our numerical scheme.

Firstly, as the imposed fluid velocity u imp depends on a provisional velocity u ⋆ (cf. Eq. ( 11)), a time splitting error is made during its calculation which may slightly affect the rate of convergence. Secondly, in the case of fluid domain Ω f totally embedded in the computational domain Ω, the calculation of the new pressure P n+1 by solving a Poisson like problem (cf. Eq. ( 21)), leads to numerical difficulties (ill-posed problem) as there is no pressure Dirichlet boundary condition imposed on the fluid domain boundary. In order to encompass this difficulty, a partial solution may consist in prescribing the pressure on one cell in the fluid domain Ω f . In doing so, we get a smooth pressure field but the velocity field is locally disturbed around the place where the pressure is prescribed.

Finally, as observed in the case of a Poiseuille flow, the graphs presented in Fig. 14 show that the calculations performed on coarse grids with the linear interpolation scheme give better results (in terms of the errors) than those done on finer grids with the base model.

[Figure 14 about here.]

Laminar flows around a cylinder

In this part, we are concerned with laminar flows around a circular cylinder. This problem has been the object of many experimental and numerical studies. In the latter case, problems involving steady and unsteady fluid flows past a static cylinder have received a particular attention. In contrast, to our knowledge, the case of a fluid flow around a rotating cylinder has been rarely investigated.

In this paper, we consider laminar fluid flows around static and rotating cylinders. The flow is characterized by the Reynolds number Re = U D ν where U is the oncoming velocity, D the cylinder diameter and ν the kinematic viscosity of the fluid. Our simulations are performed both in the steady laminar regime (i.e. Re ≤ 47 [START_REF] Guyon | Physical Hydrodynamics[END_REF]) with Re = 20 and the unsteady laminar regime with Re = 100. Tab. 2 gathers the test case parameters. We also introduce the dimensionless number β in order to characterize the rotation of the solid. It is defined as the ratio of the rotational velocity ω over the oncoming velocity and reads:

β = Dω 2U (36)
[Table 2 about here.]

[Figure 15 about here.]

As illustrated in Fig. 15, the computational domain Ω corresponds here to a square of length L with an immersed cylinder centered at the point of coordinates (0, 0). 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 [START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF].

In this work, the ratio of the length L over the diameter of the cylinder D is set to 60. This allows us to obtain, with a reasonable computational cost, results in good agreement with those given in the literature. Symmetry conditions are prescribed on Γ s while a null pressure is imposed on Γ out . As far as the boundary condition on Γ in is concerned, we consider an uniform normal velocity u in = U and a null tangential velocity. Regarding the IBCs, the following Dirichlet boundary condition is considered:

u s = ωD 2 e r on Σ (37)
where e r is the radial unit vector. Note that in the case of non-rotating cylinder (i.e. ω = 0), Eq. ( 37) reduces to a no-slip boundary condition.

Steady laminar case : Re = 20

The simulations presented and discussed herein have been performed with the linear interpolation model and the base one on uniform Cartesian grids. Qualitatively, these two models provide similar results.

Fig. 16(a) and Fig. 17(a) show respectively the streamlines and the vorticity contours obtained with a static cylinder. In this case, the flow pattern is characterized by a pair of two steady symmetric vortices attached to the surface of the cylinder. These results are in very good agreement with those proposed in the literature (e.g. [START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF][START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF][START_REF] Taira | The immersed boundary method: a projection approach[END_REF][START_REF] Linnick | A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains[END_REF][START_REF] Fornberg | A numerical study of steady viscous flow past a circular cylinder[END_REF][START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds number[END_REF]). The rotation of the cylinder disturbs the flow pattern. Indeed, as illustrated in Fig. 16(b) with the streamlines and in Fig. 17(b) with the vorticity contours, the flow becomes asymmetric with especially, a complete disappearing of the two vortices located behind the cylinder in the static case. These results are in very good agreement with [START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF] and also with those published in [START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF]. At this stage, we refer the reader to the work of Stojković & al. [40] in order to understand the effect of the rotation on the flow around the cylinder. In our case, namely for Reynolds numbers lower than the critical one (Re ≤ 47), Stojković & al. [START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF] explain that the upper vortex disappears and the lower one goes away from the surface of the cylinder with a significant decrease in size. But, in the considered regime (i.e. Re = 20)

the steady vortex obtained in the static case is small (cf. Fig. 16(a)) and thus, due to the rotation of the cylinder, it completely disappears.

[Figure 16 about here.]

[Figure 17 about here.]

To conclude the flow description, a comparison of the results presented in Fig. 18(a) and Fig. 18(b) clearly shows that the rotation also impacts the pressure distribution around the cylinder. Indeed, the pressure force rotates in the opposite direction to the cylinder which is also in good agreement with [START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF].

[Figure 18 about here.]

The results provided by the linear interpolation model and the base one are also compared quantitatively in terms of drag coefficient C d = F x /0.5U 2 D (static and rotating cases), lift coefficient C l = F y /0.5U 2 D (rotating case), recirculation length L w (static case) and direction of the total force θ = tan -1 (C l /C d ) (rotating case). This comparison is done using three levels of refinement ranging from D/h = 10 to D/h = 40 with h the size of Cartesian cells. Here-above, F x and F y correspond respectively to the tangential and normal components of the total force F t defined by:

F t = (-pI + ν∇u) • n dS ( 38
)
where n is the outward normal unit along the surface Γ.

The physical coefficients obtained with a static cylinder are summarized in Tab. 3 and are confronted with the set of data proposed in [START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF][START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF][START_REF] Taira | The immersed boundary method: a projection approach[END_REF][START_REF] Linnick | A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains[END_REF][START_REF] Fornberg | A numerical study of steady viscous flow past a circular cylinder[END_REF][START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds number[END_REF]. Fig. 19 illustrates the time convergence of these hydrodynamic coefficients.

[Table 3 about here.]

[Figure 19 about here.]

It is interesting to note that all the couples of coefficients (C d , L w /D) corresponding to the linear interpolation model are in good agreement with those presented in the literature whatever the number of cells D/h is, whereas the base model requires the finest resolution.

[Table 4 about here.]

Tab. 4 summarizes the values of C d , C l and θ obtained with the linear interpolation scheme and the base model in the rotating case. These coefficients are compared to those given in [START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF][START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF][START_REF] Ingham | A numerical investigation into steady flow past a rotating circular cylinder at low and intermediate Reynolds numbers[END_REF][START_REF] Badr | Steady and unsteady flow past a rotating cylinder at low Reynolds numbers[END_REF]. The trend observed in Tab. 4 is similar to the one previously highlighted in the static case and, as expected, the linear interpolation model yields better results than the base model. Now, a grid convergence study is conducted in order to estimate the numerical rate of convergence of the proposed method. For this purpose, without loss of generality, a smaller

computational domain is chosen, namely Ω = [-10D, 10D] × [-10D, 10D].
Resizing the computational domain allows us to consider fine grids with a reasonable computational cost. At this stage, this approach is justified because our goal is not to calculate physical coefficients but to study the behavior of our penalized direct forcing method. As no analytical solution exists, we compute the error norms ε 2 and ε ∞ (cf. Eq. ( 31)) assuming that the solution obtained on the finest mesh is the reference solution. Five levels of refinement are used ranging from h/D = 10 -1 to h/D = 6.25 × 10 -3 . For the sake of clarity, we begin by discussing the case of the static cylinder. Fig. 20 and Fig. 21 present, respectively, the evolution of the error norms ε 2 and ε ∞ measured over the whole fluid domain. As expected, Fig. 20 shows that the numerical rate of convergence in L 2 norm of the PDF method employed with the linear interpolation scheme is consistent with the second order of accuracy. The one calculated with the base model is slightly higher than one. However, as illustrated in Fig. 21, the convergence in L ∞ norm is much lower with a numerical rate of convergence close to one.

[Figure 20 about here.]

[Figure 21 about here.]

[Figure 22 about here.] Fig. 22 presents the evolution of ε ∞ calculated over 90% of the fluid domain. In this case, it is interesting to notice that the numerical rate of convergence tends toward the second order of accuracy. This indicates that the maximum of the error is located in the immediate vicinity of the immersed interface. Such a behavior has already been pointed out in the frame of the Cartesian methods by Cheny and Botella [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] which explain it by the piecewise approximation of the pressure in the Cartesian cells near the immersed interface. Actually, it is difficult for us to explain why the proposed PDF method behaves like that. Such a study will be the object of future improvements. In Fig. 20 and Fig. 22, it is also interesting to remark that the values of ε 2 and ε ∞ supplied by the linear interpolation scheme on coarse grids are lower than those measured with the base model on the finest grid. This trend is similar to the one observed in Section 5.1 and Section 5.2. Now, we focus on the results obtained with a rotating cylinder. The error norms ε 2 and ε ∞ measured over the whole fluid domain are depicted on Fig. 23 and Fig. 24, respectively. In Fig. 23, we observe that the rate of convergence in L 2 norm of the proposed method used with the linear interpolation model is close to two while it is slightly superior to one with the base model. In Fig. 24 and Fig. 25, the graphs clearly show that a calculation of the L ∞ norm over 90% of the fluid domain yields a better rate of convergence. Again, this means that the maximum of the error is very close to the immersed interface. These results confirm the trend observed with a static cylinder. Here again, the values of ε 2 and ε ∞ calculated on coarse grids with the linear interpolation scheme are lower than those computed with the base model.

[Figure 23 [2, 13-15, 38, 40, 41] in the static case and [START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF][START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF][START_REF] Kang | Laminar flow past a rotating circular cylinder[END_REF] in the rotating case).

[Figure 26 [Table 5

3-D flow induced by an ellipsoidal stirrer

In order to illustrate the ability of the PDF method with the linear interpolation model to deal with 3-D moving geometries, we describe here the flow induced by a stirrer in a cylindrical chemical reactor. The context is the vitrification process for the storage of radioactive wastes involving viscous flows (molten glass at high temperature incorporating the ultimate waste). In industrial (confidential) simulations both the vessel structure, the apparatus of measurement and the mechanical stirrer are modeled by IBCs. Here, we present a more academic test case for which the industrial stirrer is replaced by an ellipsoid in rotation around an inclined axis in the y-z plan in a cylindrical vessel. This computation was run using a pretty rough 15x15x10 Cartesian mesh. Fig. 30 shows the velocity and pressure distributions. The dynamic viscosity of the fluid is about several units. The angular velocity of the stirrer is about 60 rounds per minute, starting impulsively in a fluid at rest. Homogeneous Dirichlet BCs are prescribed on the immersed boundary describing the vessel and on the boundary of the computational domain, except on the upper surface. On this latter one, a symmetry condition is imposed. As we face to moving IBCs, an important issue concerns the treatment of the so-called "freshly fluid or solid cells". Without going into details (this will be the subject of a future paper), this is achieved by a local averaging process in the fresh fluid cells left by the moving solid. Under a CFL limitation, we set the new value of a fresh fluid cell using an arithmetic average of its neighbor cells that are in the fluid at the current and the previous time steps.

[Figure 30 about here.] For real world applications, more complex geometries are imported from Computer-Aided Design files. For instance, these computations can be useful to improve the design of the stirrer's geometry. Actually, the Trio U software has also the possibility to describe the stirrer by an IBC using in the same time a body-fitted space discretization based on the vessel.

But the use of IBCs to model this vessel allows us to drastically reduce the computational CPU time because fast solvers can be used (as multigrid solvers on Cartesian grids).

Conclusion

In this paper, we have developed a Penalized Direct Forcing (PDF) method for unsteady laminar flows around complex time-varying geometries. Such a method allows us to solve on a Cartesian grid the Navier-Stokes equations with a penalized forcing term that mimics the presence of solids through immersed Dirichlet boundary conditions. One of the purposes of this work was to develop a robust linear interpolation scheme in order to accurately reconstruct the velocity field near the immersed interfaces. The originality of our approach is that the interpolation process does not depend on a preferred direction, as often done in the literature, but is based on an averaged reconstruction of the velocity gradient around the interface and on a minimization problem that relies on a local reconstruction of the immersed interface. Regarding to the numerical scheme, the proposed method is based on a finite volume approximation with a staggered grid arrangement of the variables. Furthermore, in order to solve the governing equations, we have used a fractional step scheme that is modified in such a way that the Dirichlet immersed boundary condition for velocity is verified not only in the prediction equation but also in the correction one. Moreover, the associated homogeneous Neumann boundary conditions are also enforced on the immersed boundaries.

It is worth to notice that Neumann immersed boundary conditions could be also considered.

Several numerical experiments have been carried out using this new method. All the results demonstrate the efficiency and the potentialities of applications of our PDF method.

On the one hand, the validity and the accuracy in space of our Hybrid Cartesian/Immersed Boundary technique have been assessed throughout two academical problems involving uniform and analytic Dirichlet IBCs. From these two tests, it results that the numerical rate of convergence is (quasi-)quadratic in L 2 -and L ∞ -norms. On the other hand, the problem of a laminar flow around static and rotating cylinders has been studied both in steady and unsteady regimes. Whatever the case considered, our results are in good agreement with those published in the literature. Furthermore, a grid convergence study performed in the steady case confirms that the numerical rate of convergence of our method tends to the second order of accuracy in space.

Extension of our method to problems involving 3-D geometries is trivial and requires no coding efforts as illustrated in this paper. In fact, in the nuclear safety context, 3-D purely hydraulic numerical simulations are in progress to study vitrification processes for the storage of radioactive wastes. In addition, our method has been only assessed on 2-D problems involving static and rotating solids. In the future, it is planned to improve this assessment on 3-D problems with moving geometries. This work is already underway for the flow induced by a mixer during the vitrification process. In future works, it will be also interesting to investigate the impact of the approximation of pressure on the results and, more particularly, on the convergence of the method in the immediate vicinity of the immersed interface. Equipping our method with adaptive techniques of mesh refinement would be also useful for this task. In an other direction, we plan to work in the future on immersed wall laws for turbulent flow applications as it is often the case in the industrial context. This enhancement, that is an alternative to the mesh refinement techniques, was already mentioned by Iaccarino and Verzicco in [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] where LES computations were done using immersed boundaries and mesh enrichment techniques.

⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ ⋄ x d(x) < 0 Ω f ⋄ ⋄ ⋄ ⋄ ⋄ h Ω Ω s Ω f Σ h Figure 3: Schematic representation of solid nodes x where u imp (x) = u s (x). × x i,j v s (x i,j+1/2 ) v s (x i,j-1/2 ) u s (x i-1/2,j ) u s (x i+1/2,j ) h Ω Ω s Ω f Σ h
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × x 1 x 2 x d(x 1 ) d(x) d(x 2 ) × × × × • • • • • • h Ω Ω s Ω f Σ h
× × × × × × × × × • • • • • • • • • x Ω Ω f Σ h × × × × × × × × × • • • • • • • • • x Ω Ω f Σ h first level -1 neighboring cell cut by Σ h -1 Lagrangian facet × × × × × × × × × • • • • • • • • • x Ω Ω f Σ h x : a interfacial Cartesian node second level -2 neighbors of neighbors intersected by Σ h -2 Lagrangian facets
x associated with 3 Lagrangian facets Figure 6: Schematic representation of the collecting procedure of Lagrangian facets in the case of a Cartesian node x located in the interfacial region. In the example, three Lagrangian facets are collected. One of them is obtained during the first level where only one neighboring cells is intersected by Σ h . The last two are collected in the two "neighbors of neighbors" identified during the second level.

× × × × × ⊗ ⊗ × × • • • • • • • • • x Ω Ω f Σ h × × × × × ⊗ ⊗ × × • • • • • • • • • x Ω Ω f Σ h first level -2 neighboring cells purely fluid -2 interfacial nodes ⊗ × × × × × × × × × • • • • • • • • • x Ω Ω f Σ h x : a fluid Cartesian node second level -2 neighboring cells cut by Σ h -2 Lagrangian facets x associated with Γ i n Γ o u t Γ s Γ s Ω f Ω s Ω s Σ Σ X Y x y θ

( 9 )

 9 converges toward the Navier-Stokes body-fitted solution in L 2 (Ω) norm with O η 3/4 in the solid and O η 1/4 in the fluid (but numerical tests show O (η) whatever the subdomain is). In this work, η is set to 10 -12 .

  the fluid contribution is defined as the intersection of the nearest fluid cell with the line passing through the forcing node and the solid contribution. Another possibility, proposed by Liao & al. [12], consists in defining the fluid contribution by the nearest fluid node along a grid line direction or a diagonal direction. Similarly, the solid contribution is defined by the intersection of Σ h with the line passing through the forcing node and the fluid contribution. In the aforementioned examples, the determination of the solid contribution depends essentially on the preferred direction while the fluid contribution does not take account of the other fluid nodes in the immediate vicinity of the forcing node.

5 . 3 . 2 .Fig. 26 ,

 53226 Fig.[START_REF] Taira | The immersed boundary method: a projection approach[END_REF], Fig.27and Fig.28present, respectively, an instantaneous of the streamlines, the vorticity contours and the pressure distribution. In these figures, we remark that the flow pattern is characterized by the well-known Von Kármán vortex street which corresponds to a periodic shedding of vortex behind the cylinder. Such a phenomenon, well captured by our immersed boundary technique, is in good agreement with the results provided in the literature (e.g.[2, 13-15, 38, 40, 41] in the static case and[START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF][START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF][START_REF] Kang | Laminar flow past a rotating circular cylinder[END_REF] in the rotating case).

6 (

 6 about here.] [Figure 27 about here.] [Figure 28 about here.] Moreover, as previously done, our results are also compared in terms of hydrodynamic coefficients and confronted to those of the references formerly cited. We have computed the time averaging drag coefficient C d , the time averaging lift coefficient C l , the amplitude of their fluctuations C ′ d and C ′ l and the Strouhal number St = f D/U which is the dimensionless number used to characterize the shedding frequency f . In this work, f is estimated from the periodic variation of the lift coefficient C l . As shown in Tab. 5 (static case) and Tab. rotating case), the values of the hydrodynamic coefficients obtained with the proposed penalized direct forcing method are in good agreement with those published in the literature.

Fig. 29

 29 Fig. 29 illustrates the time evolution of these hydrodynamic coefficients.

Figure 2 :

 2 Figure 2: Schematic representation of a fictitious computational domain Ω on a Cartesian grid with the solid domain Ω s , the fluid domain Ω f and the immersed boundary Σ. Here, u s denotes the velocity of the obstacle Ω s .

Figure 4 :

 4 Figure 4: Schematic representation of the base model on a staggered grid

Figure 5 :

 5 Figure 5: Schematic representation of the linear interpolation model on a staggered grid

Figure 8 :

 8 Figure 8: Poiseuille flow in an inclined channel. Computational domain Ω with the fluid domain Ω f , the solid domain Ω s and the immersed boundaries Σ.

Figure 9 :

 9 Figure 9: Poiseuille flow in an inclined channel. Example of stationary velocity and pressure fields (right) plotted on Cartesian grid with immersed boundaries (left).

Figure 10 :Figure 12 :

 1012 Figure 10: Poiseuille flow in an inclined channel. L 2 (Ω f ) norm (left) and L ∞ (Ω f ) norm (right) of the error vs. the computational grid size h. Comparison of the linear interpolation model with the base model.

Figure 13 :

 13 Figure 13: Taylor-Couette problem. Example of the stationary velocity field (right) plotted on a Cartesian grid with immersed boundaries (left).

Figure 14 :Figure 15 :

 1415 Figure 14: Taylor-Couette problem. L 2 (Ω f ) norm (left) and L ∞ (Ω f ) norm (right) of the error vs. the ratio of the computational grid size h over the smallest radius. Comparison of the linear interpolation model with the base model.

Figure 17 :

 17 Figure 17: Vorticity contours around a rotating cylinder obtained with the linear interpolation scheme : Re = 20

Figure 18 : 20 58Figure 19 :

 182019 Figure 18: Pressure distribution around a rotating cylinder obtained with the linear interpolation scheme : Re = 20

Figure 20 :

 20 Figure 20: Flow around a static cylinder: L 2 (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over the whole fluid domain. Comparison of the linear interpolation model with the base model.

Figure 21 :

 21 Figure 21: Flow around a static cylinder: L ∞ (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over the whole fluid domain. Comparison of the linear interpolation model with the base model.

Figure 22 :

 22 Figure 22: Flow around a static cylinder: L ∞ (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over 90% of the fluid domain. Comparison of the linear interpolation model with the base model.

Figure 23 :

 23 Figure 23: Flow around a rotating cylinder: L 2 (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over the whole fluid domain. Comparison of the linear interpolation model with the base model.

Figure 24 :

 24 Figure 24: Flow around a rotating cylinder: L ∞ (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over the whole fluid domain. Comparison of the linear interpolation model with the base model.

Figure 25 :

 25 Figure 25: Flow around a rotating cylinder: L ∞ (Ω f ) norm error of the streamwise u (left) and the spanwise v (right) velocity components vs. the ratio of the computational grid size h over the smallest radius. Errors calculated over 90% of the fluid domain. Comparison of the linear interpolation model with the base model.

ω 1

 1 (s -1 ) ω 2 (s -1 ) r 1 (m) r 2 (m) ρ (kg.m -3 ) µ (Pa.s)

Table 1 :

 1 Numerical parameters used in the Taylor-Couette problem ω (s -1 ) U (m.s -1 ) D (m) L (m) ρ (kg.m -3 ) µ (Pa.s)

	Re	β

Table 2 :

 2 Numerical parameters used for the test of a laminar flow around a circular cylinder

Neumann IBCs can be also considered by interpolations involving the prescribed flux at the boundary and the velocities of the surrounding flow.

But it should be with the pressure at the previous time step or an extrapolation of the pressure as well.

This equation is similar to the velocity prediction equation involved in the first step of a non-incremental projection method.

d(x) < 0 for the solid domain and d(x) > 0 for the fluid region.

∇ h • u n+1 = 0 in Ω f and ∇ h • u n+1 = ∇ h • u n+1 imp in Ω s .We recall that for translations or solid rotations, ∇ • u n+1 imp =0.

Lagrangian facetsFigure 7: Schematic representation of the collecting procedure of Lagrangian facets in the case of a fluid Cartesian node x. In the example, two interfacial points (red ⊗) are detected in the first level of the collecting procedure and the second level results in two cells with two associated Lagrangian facets.

exact Figure 11: Poiseuille flow in an inclined channel. Plot on X = 2.9 of the ratio u(Y ) over the maximum velocity U ∞ . Comparison of the linear interpolation model with the base model in the coarse and fine grid cases.

(a) Static cylinder (b) Rotating cylinderFigure 16: Streamlines around a rotating cylinder obtained with the linear interpolation scheme : Re = 20

(a) Static cylinder (b) Rotating cylinderFigure 26: Streamlines around static and rotating cylinders : Re = 100

(a) Static cylinder (b) Rotating cylinderFigure 28: Pressure distribution around static and rotating cylinders : Re = 100

Acknowledgement

This work was granted access to the HPC resources of CINES under the allocation 2012-c2012026025 made by GENCI (Grand Equipement National de Calcul Intensif ).

List of Figures

Staggered arrangement of the unknowns on a 2-D Cartesian cell (i, j) with the CVs of pressure (blue ) and velocity components (red ). . . . . . . . 41 2 Schematic representation of a fictitious computational domain Ω on a Cartesian grid with the solid domain Ω s , the fluid domain Ω f and the immersed boundary Σ. Here, u s denotes the velocity of the obstacle Ω s . . . . . . . . . 42 3 Schematic representation of solid nodes x where u imp (x) = u s (x). . . . . . . 43 4 Schematic representation of the base model on a staggered grid . . . . . . . 44 5 Schematic representation of the linear interpolation model on a staggered grid 45 6 Schematic representation of the collecting procedure of Lagrangian facets in the case of a Cartesian node x located in the interfacial region.