C Introïni

M Belliard
email: michel.belliard@cea.fr

C Fournier

A SECOND ORDER PENALIZED DIRECT FORCING FOR HYBRID CARTESIAN/IMMERSED BOUNDARY FLOW SIMULATIONS

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

Flows around complex stationary/moving solids take an important place in lifescience context or in many engineering applications. Usually, these problems are solved by body-fitted approaches on unstructured meshes with boundary conditions directly imposed on the domain boundary. Another way is using immersed boundary (IB) techniques: the physical domain is immersed in a fixed fictitious one of simpler geometry on Cartesian grids. It allows to use efficient, fast and accurate numerical methods avoiding the tedious task of re-meshing in case of time varying geometry. In contrast, one needs specific methods to take into account the IB conditions (IBC). Here, we propose a second order penalized direct forcing method for unsteady incompressible flows with Dirichlet's IBC. It consists in adding a penalized forcing term to the initial problem, applied only on Cartesian nodes near the IB, in order to bring back the variable to the imposed one. Regarding Navier-Stokes solvers using a projection scheme, the forcing term is distributed both in the velocity prediction and in the correction equations. It leads to a natural way to prescribe the pressure boundary conditions around obstacles. Numerical experiments, performed for laminar flows around static/moving solids, assess the validity and illustrate the ability of our method, showing in particular a quadratic convergence rate.

INTRODUCTION

Fluid flows with heat and mass transfers around/inside complex stationary or moving geometries appear in a large number of situations of practical interest including biological fluid mechanic applications (e.g. blood flow in heart) or engineering applications (e.g. heat exchangers, aerospace vehicles or nuclear safety). 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 coincide with those of the physical domain. Thereby, boundary conditions are directly imposed on the physical domain boundary. However, the main drawbacks of the body-fitted like techniques lie in their lack of ability to handle complex industrial problems involving moving bodies which require the development of complex numerical schemes to deal with the difficult issue of regridding. 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 Cartesian grids. Such techniques allow to use efficient, fast and accurate numerical methods avoiding the tedious task of re-meshing caused by time-varying geometries. In contrast, as the physical boundaries are described by a set of Lagrangian points that do not generally coincide with those of the Eulerian grid or by a level-set function, numerical methods are needed to account for the Immersed Boundary conditions (IBC). Roughly, the non-boundary conforming techniques proposed in literature may be classified in 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 coeffcients 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], mimicks the presence of embedded geometries by modifying the numerical scheme in the immediate vicinity of the immersed interface. Such an approach leads to a sharp representation of the immersed interface. In the second category, rather than changing locally the numerical scheme, a forcing term is added in the governing equations. Since the Peskin's pioneering work of the Immersed Boundary Method (IBM) [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF], several IB like methods with different forcing terms (or forcing strategies) have been proposed in literature: e.g. the Goldstein's Feedback Forcing method [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF] or the Direct Forcing (DF) method 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. [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF]. The DF technique consists in directly applying the desired boundary conditions on the Cartesian nodes closest to the interface leading to a sharper representation of the interface than in the Peskin's method. In that sense, by using the terminology employed 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. 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 have gained in popularity and have 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]). More recently, Belliard & Fournier [START_REF] Belliard | Penalized direct forcing and projection schemes for Navier Stokes[END_REF][START_REF] Introni | A Second Order Penalized Direct Forcing for Hybrid Cartesian/Immersed Boundary Flow Simulations[END_REF] 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. [START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]). Links can be found with the works of Sarthou & al. [START_REF] Sarthou | Consistant velocity-pressure coupling for second-order L 2 -penalty and direct-forcing methods[END_REF] and those of Bergmann & Iollo [START_REF] Bergmann | Iollo: Modeling and simulation of fish-like swimming[END_REF].

In the present paper, after introducing the discretization of the governing equations in the Section 2, our PDF algorithm including the velocity reconstruction near the immersed boundaries is detailed in the Section 3. We have developed an original robust interpolation scheme, second-order accurate in space, that relies mainly on an averaged reconstruction of the velocity gradient near the IB and on an approximate projection operator onto the IB. Finally, in Section 4, some 2D numerical experiments are performed for laminar flows around/between solids to assess the validity and the ability of the proposed method both for IBC. Without loss of generality, we restrict our presentation to Dirichlet's IBCs. We show in particular that the numerical rate of convergence is quasi-quadratic for all studied cases.

Governing equations and numerical scheme

In this section, we focus on the governing equations and we present the resolution method, based on a fractional step method, and the space and time discretizations.

Governing equations

The governing equations used to described unsteady incompressible flows around complex obstacles Ω are the incompressible Navier-Stokes equations:

∂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 u denotes the solenoidal velocity and ν the kinematic viscosity. Here above, P is the total pressure defined by ρ∇P = ∇p -ρg where p is the hydrodynamic pressure, ρ the constant density and g the gravity force. For clarity reason, we have assumed full Dirichlet boundary conditions and we will consider f = 0 in the rest of this paper.

Space and time discretization

The time advancement of the velocity is performed by means of a degenerate fourthorder explicit Runge-Kutta scheme [START_REF] Williamson | Low-storage Runge-Kutta schemes[END_REF]. Given 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 + Σ 3 k=1 β k q k and q k = ∆tf u n + Σ k-1 m=1;k>1 β m q m + α k q k-1 with α k equals to 0, -1 2 and -2 and β k equals to 1 2 , 1 and 1 6 for k = 1, 2 or 3 respectively. At each stage k of the Runge-Kutta scheme, we get f (u n + ...) by solving the couple (u, P ) by a fractional step method, c.f. Section 2.2.1, for which the temporal discretization is based on a semi-implicit scheme: explicit discretization of the convection term and implicit discretization of the diffusion. The space discretisation is based on a finite volume approximation with a staggered grid arrangement of the primitive variables (u, P ) (velocity components at the middle of the cell edges and pressure at the cell center). In this frame, the governing equations Eqs. (1) are integrated over each control volumes ensuring the conservation of mass and momentum balance. The convection and diffusion terms are respectively approached by the QUICK and the centered schemes. We denote by ⋆ h a discrete space operator ⋆. At each stage k = 0, .., 3 of the Runge-Kutta scheme, the discrete form of the governing equations given by Eqs. (1) reads

f (u k-1 ) = u * -u k-1 ∆t in Ω (2a) u * -u k-1 ∆t + ∇ h • u k-1 ⊗ u k-1 + ∇ h P * = ν∇ 2 h u * in Ω (2b) ∇ h • u * = 0 in Ω (2c)
with suitable boundary on ∂Ω, u k = u n +Σ k m=1 β m q m , u k-1 = u n for k = 1 and (u n+1 , P n+1 ) = (u k , P * ) for k = 3 .

A projection algorithm

The fractional step method or projection method was introduced by Chorin and Temam in 1968 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, 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 incompresible 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 other. 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 our case, the first step of the algorithm consists in solving a predicted velocity u without pressure gradient term as follows 1 :

u -u n ∆t + ∇ h • (u n ⊗ u n ) -ν∇ 2 h u = 0 in Ω (3a) u = u n+1 = u D on ∂Ω (3b)
The second step of the algorithm 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 Ω (4a) ∇ h • u n+1 = 0 in Ω (4b) n • u = n • u n+1 = n • u D on ∂Ω (4c)
Finally, the correction Eqs. (4) allows to compute the new velocity u n+1 as follows

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

Penalized Direct Forcing

As previously mentioned, the immersed boundary technique consists in immersing the physical domain in a fixed fictitious one. Here above, we will denote by Ω f the physical domain (fluid) and by Ω the fictitious domain discretized on a Cartesian mesh. We denoted by Ω s the external domain (solid) defined by Ω s = Ω\Ω f . The domains Ω f and Ω s are separated by an immersed boundary Σ discretized by a set of Lagrangian points. In the numerical experiments presented in the present paper, the interface tracking is performed by means of a Front-Tracking technique, c.f. for instance [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]. We will use subscript h for discrete quantities. To begin, we introduce in a first part the concept of the penalized direct forcing method and we present a fractional step method consistent with the penalized direct forcing [START_REF] Belliard | Penalized direct forcing and projection schemes for Navier Stokes[END_REF]. Then, in a second part, we focus on the interpolation schemes in order to reconstruct accurately the velocity field (second order in space) near the immersed boundary. A detailed presentation can be found in [START_REF] Introni | A Second Order Penalized Direct Forcing for Hybrid Cartesian/Immersed Boundary Flow Simulations[END_REF].

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 governing equations used to described unsteady incompressible flows around complex geometries are given by:

∂u ∂t + ∇ • (u ⊗ u) + ∇P -ν∇ 2 u = F in Ω, (6a)
∇ • u = 0 in Ω. (6b)
Eqs. ( 6) are the incompressible Navier-Stokes Eqs.

(1) with a penalized forcing term F here defined by [START_REF] Belliard | Penalized direct forcing and projection schemes for Navier Stokes[END_REF]:

F = αχ s η (u i -u) with α > 0, 0 < η ≪ 1 and χ s → [0, 1] (7)
where χ s is the characteristic function of the obstacle, η a penalty coefficient (in this work η = 10 -12 ) and u i the imposed fluid velocity around/inside the obstacle Ω s depending on the obstacle velocity u s but also on the fluid velocity u following the requested accuracy. Eq. ( 7) can be view as an implicit limit version of direct forcing expression [START_REF] Mohd-Yusof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF][START_REF] Introni | A Second Order Penalized Direct Forcing for Hybrid Cartesian/Immersed Boundary Flow Simulations[END_REF]. The forcing term F is only applied on Cartesian nodes near or inside the immersed boundary (i.e. χ s (x) > 0), leading to u i (x) = u(x). Outside, the classical incompressible Navier-Stokes equations are recovered in the fluid domain Ω f (i.e. χ s (x) = 0).

Resolution of the Navier-Stokes equations: a consistent fractional step method

In this section, we propose a fractional step method to solve the Navier-Stokes equations Eqs. [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF] in the framework of the Penalized Direct Forcing. The new feature of our algorithm is that the forcing term is distributed both in the prediction and the correction stages of the projection. This leads to a natural consistent scheme in the sense that the immersed boundary conditions are well satisfied not only by the predicted velocity in the prediction stage but also by the new velocity at the end of the projection stage [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF].

As usual, the first step of our scheme consists in solving a predicted velocity u, see Eqs.

(3), with the forcing term given by Eq. ( 7) and α = 1 ∆t :

u -u n ∆t + ∇ h • (u n ⊗ u n ) -ν∇ 2 h u = χ s η∆t (u n+1 i -u) in Ω (8)
In this equation, the imposed fluid velocity u n+1 i depends on the obstacle velocity and on the fluid velocity around the obstacle. This dependency is here treated in an explicit timediscretization way, solving first a free-obstacle explicit N.-S. equations ( u ⋆ ). Then the imposed velocity u n+1 i is expressed in terms of u n+1 s and u ⋆ . The solution of Eq. ( 8) respects the Dirichlet boundary conditions on Σ h up to a given order depending on the space interpolation scheme, see Section 3.3. The second step consists in a correction stage that reads as usual, see System (4), except for the Eq. (4a):

u n+1 -u ∆t = -∇ h P n+1 + χ s η∆t ( u -u n+1 ) in Ω (9)
that can be simplified as

ρu n+1 -u ∆t = -∇ h P n+1 in Ω. ( 10
)
This equation is similar to (4a) with a modified density ρ = (1 + χ s /η). Therefore the rest of the projection algorithm is similar to the standard one using this modified density ρ. In fact, we get a consistent fractional step scheme in the sense that the Dirichlet immersed boundary condition is well satisfy by the new velocity u n+1 (not only by u) and that the homogeneous Neumann immersed boundary condition is satisfy by the pressure on the obstacles.

Reconstruction of 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 into account the immersed boundary conditions. An simple 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] (referred here as the base model, first-order accurate in space). This model leads to a stepwise description of the immersed interface Σ h .

Among the numerical algorithms proposed in literature to reconstruct the velocity field in the interfacial region, the most widely employed to improve the accuracy of the solution in the neighborhood of Σ h is based on an interpolation or an extrapolation procedure (e.g. [START_REF] Peller | High-order stable interpolations for immersed boundary methods[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[END_REF][START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF][START_REF] Liao | Simulating flows with moving rigid boundary using immersed-boundary method[END_REF][START_REF] Gilmanov | A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies[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]) involving solid and fluid velocities to calculate the velocity at the forcing nodes. Unlike the classical approaches, consisting in a local interpolation or extrapolation along a specific direction (often the grid-line directions or the orthogonal direction onto Σ h ), our proposed linear interpolation method, second-order accurate in space, involves an averaged reconstruction of the velocity gradient near the immersed boundary [START_REF] Introni | A Second Order Penalized Direct Forcing for Hybrid Cartesian/Immersed Boundary Flow Simulations[END_REF]. The fluid contribution is built following similar ideas used in [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF] and the solid contribution is determined by means of a minimization problem. It consists in estimating the imposed fluid velocity u i at the forcing or penalized node x according to

u i (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 (11)
where x p is the p-th fluid node in the immediate vicinity of the penalized node x. Hence, we take into account the local influence of the fluid flow around x. Π Σ (x) denotes the projection of x onto the immediate neighborhood of Σ h and is defined throughout an algorithm based on the following minimization problem:

Find z (= Π Σ (x)) ∈ Ω such that J(z) = inf y∈Ω J(y) (12)
where J(•) is defined by ∀ y ∈ Ω, J(y) = yx 2 . To tackle this minimization problem, we partially reconstruct the immersed boundary in the immediate vicinity of x, collecting all the Lagrangian facets m that belong to the cells intersected by Σ h around x and determining for each of them his plan equation d j=1 c mj y j = f m . Here d is the spatial dimension of the problem, y j the coordinates of one vertex belonging to the m-th facet and c mj the coordinates of the unit normal vector at the centroid of this facet. Then for all the facets associated with the Cartesian node x, we have C • y ≈ f. To solve this minimization problem, we use an Uzawa algorithm which can be summarized as follows:

1. Initialisation (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
and λ k the vectorial Lagrange multiplier at k-th iteration associated with the constraint C • y k+1 = f. 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, leading to a prohibitive CPU-time cost. To overcome this difficulty, we only selects the most representative facets and detects also those that are collinear in order to remove the duplicate. This algorithm is very robust and z = Π Σ (x) ∈ Ω converges toward z ∈ Σ h when h is decreasing.

Numerical experiments

This section is devoted to the numerical validation of the proposed Penalized Direct Forcing method. On the one hand, Poiseuille flow in an inclined channel, cylindrical Couette flow and steady flow around static/rotating cylinder are considered to assess the validity and the ability of our method both for Dirichlet's IBC. Grid convergence studies have been done in order to obtain the numerical rate of convergence of the method in the L

2 (Ω f ) norm ε 2 and the L ∞ (Ω f ) norm ε ∞ given by ε 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 |
where the upperscript i denotes the i-th face of the Eulerian grid, u i ref a referenced velocity calculated on the i-th face and N f the total number of faces in the fluid region Ω f . All these academical numerical tests have been performed with both the linear interpolation scheme and the base model. On the other hand, interestingly for pratical purposes, 2D calculations involving unsteady laminar flow and the interpolation scheme are presented.

All the calculations have been performed by using the CFD code Trio U developed at the French "Commissariat à l' Énergie Atomique et aux Énergies Alternatives" [START_REF] Calvin | A development framework for parallel CFD applications: Trio U project. Su-perComputing in Nuclear Applications (SNA-2003)[END_REF].

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). The computational domain Ω is define by Ω s ∪Ω f and the immersed boundary by Σ = Ω s ∩ Ω f . Null pressure boundary conditions are imposed on ∂Ω except on the entry where we prescribed the following parabolic velocity

u in = U ∞ -Y 2 , 0 T with Y = y-x √
2 , (X, Y ) the system of coordinates in the frame associated with the inclined channel and U ∞ = 0.605 m.s -1 the maximum velocity. For this numerical test, no-slip boundary condition are imposed on the immersed boundary Σ (i.e. u s = 0). For a convergence study, four computational grid sizes have been considered over the range 25.10 -3 ≤ h ≤ 2.10 -1 . Fig. 1 represents the evolution of the L 2 (Ω f ) norm ε 2 and the L ∞ (Ω f ) norm ε ∞ of the error. As expected, the proposed linear interpolation model leads to a quasiquadratic numerical rate of convergence while the numerical order of the method using the base model is about one. We may remark that, as the whole, the errors in the case of the linear interpolation model are at least one order of magnitude lesser than the ones obtained with the base model. Indeed, the approximation obtained in the case of the coarsest grid with the linear interpolation model is very close to the one computed with the base model on the finest grid. Therefore, for this test, the results clearly show that the proposed method coupled with the linear interpolation model allows to obtain an accurate solution at lower computational cost.

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]). The computational domain is defined by a square

Ω = Ω s ∪ Ω f = [0, L] × [0, L]
in which are immersed the boundaries Σ 1 and Σ 2 mimicking the inner (r 1 = 0.1m) and outer (r 2 = 0.2m) cylinders, respectively. The length L is about 2 (r 1 + r 2 ). The inner cylinder rotates clockwise (ω 1 = 1s -1 ) while the outer cylinder rotates in the counterclockwise (ω 2 = -1s -1 ). In overall calculations, we assume that ρ = 1kg.m -3 , µ = r 1 /ℜ Pa.s and that the Reynolds number, defined by ℜ = |ω 1 |r 2 1 /ν, is set to be 1. Such an assumption allows us to write Ta < Ta c where Ta = 3/2 is the Taylor number defined by Ta = 0.5ℜ 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]. This relation implies that the following simulations remain strictly planes. Regarding the boundary conditions, we prescribed symmetry conditions on the boundaries of the computational domain ∂Ω and the following analytical velocities on the immersed boundaries Σ 1 and Σ 2 : ∀ i = 1, 2 u s = ω i r i e r on Σ i where e r is the radial unit vector. The flow is initially at rest and converges toward a steady state. Here again, a grid convergence study have been done in order to estimate the numerical order of the proposed method. Five meshes have been considered over the range 25.10 -3 ≤ h/r 1 ≤ 10 -1 where h denotes the step size. Fig. 2 presents the L 2 (Ω f ) norm ε 2 and the L ∞ (Ω f ) norm ε ∞ of the error. As in the case of the Poiseuille flow, a quasi-linear numerical rate of convergence is obtained with the base model and a quasi-quadratic rate with the linear interpolation scheme. Finally, as previously observed, the graphs presented in Fig. 2 show that the calculations performed on the coarsest grid with the linear interpolation scheme give better results (in terms of the error magnitude) than those of the finest grid with the base model.

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 receive a particular attention. In contrast, to our knowledge, the case of fluid flows around rotating cylinder have been rarely investigated.

In this paper, we consider laminar fluid flows around static and rotating cylinders. Lets U = 1 m.s -1 be the oncoming velocity, D = 1 m the cylinder diameter and ν = 1/ℜ m 2 .s -1 the kinematic viscosity of the fluid. Our simulations are performed both in the steady laminar regime (i.e. ℜ = U D ν ≤ 47 [START_REF] Guyon | Physical Hydrodynamics[END_REF]) with ℜ = 20 and in the unsteady laminar regime with ℜ = 100 (ρ = 1kg.m -3 and µ = 1/ℜ Pa.s). We also introduce the dimensionless number β in order to characterize the rotation of solid. It is defined as the ratio of the rotational velocity ω over the oncoming velocity and reads β = Dω 2U . Here, the rotational velocity ω is equal to 2 s -1 (β = 1). The computational domain Ω corresponds to a square of length L with an immersed cylinder centered at the coordinate origin. 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 literature. Symmetry conditions are prescribed on ∂Ω except at th e in-flow boundary (uniform normal velocity u in = U and null tangential velocity) and at the out-flow one (null pressure). Regarding the immersed Dirichlet boundary conditions on Σ, we have: u s = ωD 2 e r where e r is the radial unit vector. In the case of a non-rotating cylinder (i.e. ω = 0), it boils down to a no-slip boundary condition.

Steady laminar case: ℜ = 20

Qualitatively, the linear interpolation model and the base one provide similar results concerning the flow pattern. Fig. 3(a) shows the streamlines 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 in agreement with the results proposed in literature (e.g. [START_REF] Taira | The immersed boundary method: a projection approach[END_REF][START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF][START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[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 as illustrated in Fig. 3(b). The flow becomes asymmetric with especially, a complete disappearing of the two vortices located behind the cylinder in the static case. Again, 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 concur 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. [34] in order to understand the effect of rotation on the base 0.82 0.98 0.925 0.92 0.9 0.94 0.93 0.91 linear 0.92 0.91 0.9

Table 1: Hydrodynamic coefficients associated with the problem of steady flow around static cylinder flow around the cylinder. Quantitatively, the linear interpolation model and the base one are compared 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). Here, F x and F y correspond respectively to the tangential and normal components of the total force F t . This comparison is done by using three levels of refinement ranging from D/h = 10 to D/h = 40 with h the size of Cartesian cells. The physical coefficients obtained with a static cylinder are summarized in Tab.1 and confronted with the set of data proposed in [START_REF] Taira | The immersed boundary method: a projection approach[END_REF][START_REF] Choi | An immersed boundary method for complex incompressible flows[END_REF][START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[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]. Whatever the number of cells D/h is, the linear interpolation model is in good agreement with literature concerning the evaluation of the couples of coefficients and are compared with those given in [START_REF] Stojković | Effect of high rotation rates on the laminar flow around a circular cylinder[END_REF][START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[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 is similar to the one previously noted in the static case and, as expected, the linear interpolation model yields better results than the base model.

(C d , L w /D),
Restricting ourselves to the behavior of the penalized direct forcing method, a smaller com- to perform a grid convergence study, allowing to consider fine grids with a reasonable computational cost. Finest grid solution is chosen as the reference solution to compute the error norms ε 2 and ε ∞ . Five levels of refinement are used ranging from h/D = 10 -1 to h/D = 6.25 × 10 -3 . Fig. 4 and Fig. 5 presents, respectively, the evolution of error norms ε 2 and ε ∞ measured on the whole fluid domain. As expected, the numerical rate of convergence in L 2 norm of the method employed with the linear interpolation scheme is consistent with the second order in space and the one calculated with the base model is slightly higher that one. However, the convergence in L ∞ norm is lower with a numerical rate of convergence close to one for the both approaches. But, as previously point out, the linear interpolation scheme on the coarsest grid is again more accurate than the base model on the finest grid. In fact, if we compute the ε ∞ norm on 90% of the fluid domain (far from the immersed boundary), the numerical rate of convergence tends to the second order for the linear interpolation scheme (1.87 for u and 1.84 for v). This indicates that the maximum of the error is located in the immediate vicinity of the immersed interface, as already pointed out in the frame of 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]. Concerning the case of the rotating cylinder, Fig. 6 and Fig. 7 present the evolution of the ε 2 and ε ∞ norms measured on the whole fluid domain. These results confirm the trend observed with a static cylinder. The error norm ε ∞ is about 1.87 for u et 1.84 for v.

Unsteady laminar case: ℜ = 100

Here we face the simulation of the unsteady regime with ℜ = 100, using only the linear interpolation model. The geometrical features are identical to those considered previously in the steady case and all the calculations have been done by using a non-uniform grid with fifty Cartesian cells in the diameter of the cylinder. Fig. 8 presents an instantaneous view of vorticity contours. The flow pattern is characterized

CONCLUSIONS

In this paper, we have presented the Penalized Direct Forcing method, using here a finite volume space approximation with a staggered grid arrangement of variables. Such a method allows to solve on a Cartesian grid the Navier-Stokes equations with a penalized forcing term taking into account immersed boundary conditions (here Dirichlet BC, but it is possible to deal with Neumann BC). We have developed an original robust linear interpolation scheme in order to reconstruct the velocity field near the immersed boundary, not depending on a preferred direction, as often done in literature, but based on a minimization problem relying on a local reconstruction of the immersed boundary. Furthermore, the resolution of governing equations is done by using a fractional step scheme that is modified in such a way that the immersed boundary condition for velocity is verified not only in the prediction equation but also in the correction one. Several numerical experiments have been carried out. All the results demonstrate the effectiveness and the potentiality of application of our method. First, the validity and the accuracy in space have been assessed throughout two academical problems: the numerical rate of convergence is quasi-quadratic in L 2 -and L ∞ -norms. Second, the problem of laminar flows around static and rotating cylinders has been studied both in the steady and unsteady regimes. Our results are in good agreement with those published in literature and a steady-case grid convergence study confirms the numerical rate of convergence of the method. The extension to problems involving 3D geometries is trivial and requires no coding efforts. Actually, in the nuclear safety context, 3D purely hydraulic numerical simulations are in progress to study the vitrification processes for the storage of radioactive wastes handling problems with moving geometries. In this case, an important issue concerns the treatment of the so-called "freshly fluid or solid cells".

Figure 1 :

 1 Figure 1: 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 2 :

 2 Figure 2: 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 3 :

 3 Figure 3: Streamlines around rotating cylinder obtained with the linear interpolation scheme: ℜ = 20

Figure 4 :

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

Figure 5 :

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

Figure 6 :

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

Figure 7 :

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

Figure 8 :

 8 Figure 8: Vorticity contours around static and rotating cylinders: ℜ = 100

Table 2 :

 2 whereas the finest grid resolution is required by the base model. Tab.2 summarizes the obtained values of C d , C l and θ in the rotating case • 58.97 • 57.68 • 54.31 • ∼ 56 • 53.66 • 53.87 • linear 53.95 • 55.95 • 56.05 • Hydrodynamic coefficients associated with the problem of steady flow around rotating cylinder

			D/h			References
		10	20	40	[2]	[34]	[39]	[40]
	C d	base 1.8968 1.8703 1.8608 1.888 ∼ 1.85 1.925 2.000 linear 1.9104 1.8746 1.8679
	C l	base 3.0284 3.1097 2.9419 2.629 ∼ 2.75 2.617 2.740 linear 2.6248 2.7740 2.7745
	θ	base 57.93

Table 3 :

 3 Hydrodynamic coefficients associated with the problem of unsteady flow around static cylinder

		Present						References
		results	[34]	[36]	[2]	[13]	[12]	[32]	[14]
	C d 1.347	1.3371	1.34	1.392	1.34	1.35	1.317	1.3757
	C C d ′ ′ l	±0.009 ±0.0091 ±0.009 0.326 0.3259 0.333	--	±0.011 ±0.012 ±0.009 ±0.0096 0.315 0.303 0.349 0.3393
	St 0.165	0.165	0.166 0.164	0.167	0.170	0.1692
						Present		References
						results	[34]	[2]	[41]
			C d	1.12	1.1080	1.189	1.0979
			C	′ d	±0.11 ±0.0986 ±0.1195 ±0.0988
			C l	2.51	2.504	2.405	2.4833
			C	′ l	±0.37 ±0.3616 ±0.4427 ±0.3603
			St 0.165	0.1658	0.1732	0.1650

Table 4 :

 4 Hydrodynamic coefficients associated with the problem of unsteady flow around rotating cylinder

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


M Belliard ♦,

M Chandesris ♦

J Dumas ♦

Y Gorsse ♦,

D Jamet ♦

C Josserand

ON THE CONTROL OF SPURIOUS FORCE OSCILLATIONS FOR MOVING BODY PROBLEMS USING AN IMMERSED BOUNDARY METHOD

Keywords: Immersed Boundary Method, Moving Bodies, Spurious Force Oscillations, Fluid Structure Interaction

We present a Direct Forcing (DF) approach for simulations of flows over moving bodies, on Cartesian grids, capable of reducing the amplitude of the force oscillations, by using a regularized approach.

INTRODUCTION

Fluid-structure interaction (FSI) problems are tricky to simulate due to the complexity of the physical phenomena and the presence of time varying geometries. In the well-known body-fitted approach the boundary conditions (BC), that are critical for FSI issues, are exactly imposed on the physical domain boundary. However accurate this approach may be, it is sometimes difficult to deal with the re-meshing issue, specially for large fluidmesh deformations. Another approach consists in using immersed boundary methods (IBM) [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF]. The general idea is to consider only a fluid domain in which the solid boundaries are immersed. We recover the immersed BC by adding a supplementary term in the governing equations of the fluid, referred to as the forcing term. Here, we consider the Direct Forcing (DF) method [START_REF] Mohd-Yussof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF] that is particular discrete version of the IBM for rigid solids. This method has been successfully applied to various problems [START_REF] Fadlun | Combined immersedboundary finite-difference for three-dimensional complex flow simulations[END_REF]. Unfortunately, one undesirable property is the generation of spurious force oscillations (SFO) [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] 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. Here, we present a Regularized Direct Forcing (RDF) approach [START_REF] Belliard | A regularized formulation of direct forcing IBM to control spurious force oscillations in order to simulate complex fluid structure interaction problems[END_REF] for simulations of flows over moving bodies capable of reducing the amplitude of the SFO. We have identified the primary source of these oscillations as the temporal discontinuity of the forcing term across the fluid-solid interface, concerning mainly the dead cells. Indeed, when the velocity point comes near the solid boundary, the original DF approach consists in forcing the previously computed fluid velocity to a prescribed one. Therefore the forcing term is highly discontinuous in time for dead cells leading to the famous spurious force oscillations. In the proposed RDF approach, this prescribed fluid velocity is chosen in a such a way that the magnitude of the forcing term tends to zero with the distance to the immersed interface (or the solid volume fraction). By modifying the definition of the prescribed fluid velocity, it is possible to change the order of approximation of the method. A first order and a second order version of the method are presented here. The results obtained with the present RDF approach are similar to those obtained with much more complicated method [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][START_REF] Lee | An implicit ghost-cell immersed boundary method for simultions of moving body problems with control odd spurious force oscillations[END_REF]. Moreover this method is very easy to implement, effective and improve the numerical precision by comparison with the original DF approach.

NUMERICAL METHODS

The governing equations used to describe unsteady incompressible flows are the incompressible Navier-Stokes (NS) equations -here with Dirichlet BC for the sake of the presentation-:

∂u ∂t + ∇ • (u ⊗ u) + 1 ρ ∇p -ν∇ 2 u = 0 in Ω ∇ • u = g in Ω (1) u = u D on ∂Ω
where Ω is the computational domain, ∂Ω its boundary, u the fluid velocity, ν the kinematic viscosity, p the pressure, ρ the fluid density and g bulk forces. A non-incremental fractional-step scheme is used to solve the governing equations (1) [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].

Following the IBM developed by [START_REF] Mohd-Yussof | Combined Immersed Boundaries/B-Splines Methods for Simulations of Flows in Complex Geometries[END_REF], the computational domain Ω corresponds to a uniform Eulerian grid including the fluid domain Ω F and the embedded solid domain Ω S such as Ω = Ω F ∪ Ω S . The presence of embedded time-varying geometry Ω S is taken into account thanks to a source term F DF , added to the time semi-discretized form of the Navier-Stokes Eqs. [START_REF] Peskin | Numerical Analysis of blood flow in heart[END_REF]. For instance, with an explicit Euler scheme and considering a projection solver, this is done in two steps:

u ⋆ -u n ∆t + ∇ • (u n ⊗ u n ) + ∇P n -ν∇ 2 u n = 0 in Ω, (2)
with u ⋆ a predicted velocity without taken into account the immersed boundary Σ. This step depends on the time scheme. The forcing step is:

ũ = u ⋆ + ∆tF DF ( 3
)
where ũ is the predicted velocity that satisfied the IBC and ∆t the time step. Finally, the projection and the correction steps remains classical:

1 ρ ∆p n+1 = 1 ∆t ∇ • ũ in Ω ( 4
)
u n+1 = ũ - ∆t ρ ∇p n+1 in Ω (5)
The source term F DF is defined by:

F DF = χ m u i -u ⋆ ∆t in Ω (6)
where u i is the imposed fluid velocity and χ m is the phase indicator function of the imposed velocities domain Ω m I :

χ m (x) = 1 for x ∈ Ω m I 0 otherwise. ( 7
)
where x is the position. The imposed fluid velocity can obtained by interpolations of various orders from the BC and fluid velocity. Here we consider interpolations of order 1 (base model: the imposed velocity is equal to the solid one u D ) and of order 2 (linear model: linear interpolation [START_REF] Introïni | A second order penalized direct forcing for hybrid cartesian/immersed boundary flow simulations[END_REF]).

The space discretization is based on a finite volume approximation 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. The governing Eqs. (1) are integrated over each volume control ensuring the conservation of mass and momentum balance. The convection and diffusion terms are respectively approached by the QUICK and centered schemes.

ANALYZE OF THE SFO

Let us consider the 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] as sandbox. It consists in a circular cylinder of diameter D oscillating along the x-direction in a fluid at rest:

x c (t) = x c (0) + X 0 (1 -cos(2πf 0 t)) (8) y c (t) = y c (0)
where (x c (0), y c (0)) are the initial coordinates of the cylinder center, X 0 the oscillations amplitude and their frequency f 0 , see the Fig. 1 Mainly it is the pressure drag component that is disturbed by the spurious oscillations whereas the time history of the friction drag is almost regular. Therefore, the SFO are in fact due to the spurious pressure oscillations (SPO). It is also interesting to notice that C DF is nearly equal to C T .

Analyzing the SPO origin in the case of the DF base model, we obtain under reasonable hypothesis [START_REF] Belliard | A regularized formulation of direct forcing IBM to control spurious force oscillations in order to simulate complex fluid structure interaction problems[END_REF] that the measure of the DF term variation between the time steps n and n + 1 can be approached by:

Ω n+1 s,h F n+1 DF dΩ - Ω n s,h F n DF dΩ = O(1) K∈Ω n+1 s,h ∩Ω n s,h + O(∆x) + O( ∆x 2 ∆t ) K∈DC + O(∆x) K∈F C (9) = O(1) + O(∆x) + O( ∆x 2 ∆t )
where ∆x is the space step, Ω n * s,h is imposed velocities domain at time n * , the F C (fresh cells) is the set of velocities that became freshly fluids (K ∈ Ω n s,h Ω n+1 s,h ) and DC (dead cells) the set of velocities that became freshly forced (K ∈ Ω n+1 s,h Ω n s,h ). This equation implies that the SPO decrease with decreasing the grid spacing and increasing the computational time step. The exponents are in good agreement with [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF]. Additional studies conclude that the DC are the main source of pressure oscillations [START_REF] Belliard | A regularized formulation of direct forcing IBM to control spurious force oscillations in order to simulate complex fluid structure interaction problems[END_REF].

REGULARIZED DIRECT FORCING

The classical DF method consists in modeling χ m in consistence with a sharp transition of the forcing term at the forced-free interface, cf. Eq. [START_REF] Belliard | A regularized formulation of direct forcing IBM to control spurious force oscillations in order to simulate complex fluid structure interaction problems[END_REF]. Then when a fluid point is reached by the solid boundary, we shapely jump from a velocity computed by the Navier-Stokes numerical scheme to a geometrical interpolation. Then the forcing term is highly discontinuous in time for the dead cells.

The key idea of the regularized DF formulation is to smooth this transition by replacing χ m with a smoother function. The following formulation is proposed for the imposed velocity:

u i = τ m u m + (1 -τ m )u ⋆ (10)
where τ m is a smooth function of the distance to the immersed boundary Σ. For instance, τ m can be the fraction of the forced-domain volume to the control volume associated with the fluid velocity. Hence we force the velocity proportionally to τ m in order to have a smooth transition: the fluid velocity will be forced progressively to the forced velocity as it enters in the forced domain. This formulation does not increase the computational time cost and is very easy to implement. We named it "regularized" (RG) DF because the imposed velocity is a linear combination of the NS computed velocity and the imposed interpolated velocity leading to a imposed regularized velocity.

For your tests we choose the following model for τ m :

u i =    u m if τ m = 1 τ m u m + (1 -τ m )u ⋆ if 0 < τ m < 1 u ⋆ if τ m = 0 (11)
with τ m computed thanks to the distance function. It varies linearly with the distance from 0 to 1 across a layer of thickness ∆x centered on the forced-free interface, see Fig. 2. Applying similar analysis than done in Section 3, we now obtain that the measure of the DF term variation between the time steps n and n + 1 can be approached by:

Ω n+1 s,h F n+1 RG dΩ - Ω n s,h F n RG dΩ = O(1) + O(∆t) + O(∆x) K∈Ω n+1 s,h ∩Ω n s,h + O(∆t) + O(∆x) K∈DC + O(∆t) K∈F C = O(1) + O(∆t) + O(∆x) (12)
Hence, the O( ∆x 2 ∆t ) term of Eq. ( 9) is reduced to O(∆x).

NUMERICAL APPLICATIONS

To demonstrate the effectiveness of the present DF forcing approach in comparison with the original one, we simulate two classic moving-body problems from the literature by varying both the grid spacing and the computational time step. The first test case is from [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] where we consider an oscillating cylinder in a fluid at rest. Concerning the second test case, it is an oscillating cylinder in a cross flow [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF].

The force oscillations are almost entirely suppressed for the first test case, cf Fig. 3. We used the pressure 2δ-discontinuity C 2δ P [START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF] to quantitatively describe the SPO:

C 2δ P = C n+1 P -2C n P + C n-1 P ( 13
)
where n is the time-step index. Space-step and time-step convergence studies of the amplitude of the SPO, cf. Fig. 4, conclude to a big damping of the SPO and corroborate the theoretical analysis. Also the present DF approach decreases the amplitude of the oscillations by almost one order of magnitude for the second one.

  (a).Fig. 1(b) shows the evolution over a period T 0 of the dimensionless coefficients concerning the pressure drag C P , the friction drag C ν , the total drag C T and the total contribution of the direct forcing term C DF = RRR Ω S ρF DF •e x dΩ F ref •e x with F ref a scaling vector and e x is the x-direction unit vector.

Figure 1 :

 1 Figure 1: Mittal's test case 1(a) geometry and computational grid. Pysical dimensionless coefficients time histories: C P , C ν , C T and C DF .

Figure 2 :

 2 Figure 2: Representation of χ m and τ m for the classical DF method and the Regularized DF one.

Figure 3 :

 3 Figure3: Time evolution of the pressure drag for an oscillating cylinder, test case from[START_REF] Seo | A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations[END_REF], using the original method and the present method, both with the first order version.

( a )

 a Grid spacing convergence of the oscillations (b) Time-step convergence of the oscillations

Figure 4 :

 4 Figure 4: Cp 2δ Max for the base and the linear interpolation schemes for both the DF and RG formulations. The CFl number is set to 1.

CONCLUSIONSConsequently, by choosing carefully the grid spacing and the time step it is possible to control and manage the magnitude of the force oscillations, leading possible to simulate full fluid structure interaction problems such as flow induced vibration.


Arturo Escobedo

Anne Spalanzani
email: anne.spalanzani@inria.fr

Christian Laugier
email: christian.laugier@inria.fr

Using social cues to estimate possible destinations when driving a robotic wheelchair

Keywords: Intention estimation, human aware navigation, semi-autonomous navigation, robotic wheelchair, face control, voice control

Approaching a group of humans is an important navigation task. Although many methods have been proposed to avoid interrupting groups of people engaged in a conversation, just a few works have considered the proper way of joining those groups. Research in the field of social sciences have proposed geometric models to compute the best points to join a group. In this article we propose a method to use those points as possible destinations when driving a robotic wheelchair. Those points are considered together with other possible destinations in the environment such as points of interest or typical static destinations defined by the user's habits. The intended destination is inferred using a Dynamic Bayesian Network that takes into account the contextual information of the environment and user's orders to compute the probability for each destination.

I. INTRODUCTION

Ensuring proper living conditions for an ever growing number of elderly people is a significant challenge in many countries. In this context, a very relevant application is to assist people with reduced mobility to preserve some of their independence.

Electric wheelchairs operated by a joystick are the most common mobility assistance device. However, elderly people often spend a large period of time learning to operate accurately those devices and even expert drivers loose the freedom of one of his hands in order to manipulate the joystick. This restricts the user's capacity to perform other task such as holding things, opening a door etc. Under the consideration that it is a natural human behavior to look where we are going [START_REF] Bernardin | Gaze anticipation during human locomotion[END_REF], in this paper we propose the use of the direction of the face together with a voice recognition system as input methods to drive the robotic wheelchair.

A safe movement must be guaranteed by this type of robotic wheelchairs specially because they are intended to be used by fragile persons. Other important aspects that must be considered in the design are: The system must be compliant with the user's intentions, it has to behave in a socially acceptable way and it must interact in an intuitive way with the driver.

Our approach is based on and extends the work presented in [START_REF] Escobedo Cabello | Multimodal Control of a Robotic Wheelchair: Using Contextual Information for Usability Improvement[END_REF]. The destination inference method was improved by considering non-static destinations to join groups of people in a socially acceptable way. The notion of personal space and interaction space is used by the navigation system to prevent uncomfortable situations when the wheelchair is moving among humans [START_REF] Rios-Martinez | Understanding human interaction for probabilistic autonomous navigation using Risk-RRT approach[END_REF].

We propose to estimate the user's intention in order to reduce the number of necessary commands to drive the robotic wheelchair and deal with ambiguous or inaccurate input interfaces. In this way, the wheelchair can be in charge of some part of the navigation task and alleviate the user involvement.

The problem is approached from the hypothesis that it is possible to learn typical destinations where the user spends most of his time and uses this information to predict his desired goal when driving the robotic wheelchair. Those typical destinations may be learned and then used to predict motion on the basis of a user model.

A probabilistic framework is used to model the existent relationship between the intention of the user and the observed command. The main originality of the approach lies in modeling explicitly the intentions as typical destinations and then use it to check the reliability of the user's command to decide how much preeminence it should be assigned in the final robot control.

The structure of this paper is as follows: Section II offers an overview of related works. A general description of the system architecture is presented in III; section III-A focuses on the estimation of intended destination and section III-B explains the module to detect social interactions and compute adequate meeting points to be taken as possible destinations. In section IV our experimental setup is presented and an analysis of the results is provided. Section V presents conclusions about the work and perspectives for future improvements of the system.

II. RELATED WORK

A. Robotic Wheelchair Navigation Systems

Research in robotic wheelchairs can be classified in three different categories; robotic wheelchairs with fully autonomous navigation, semi-autonomous wheelchairs and manually controlled wheelchairs with collision avoidance.

Robotic wheelchairs with fully autonomous navigation are normally designed as mobile robots equipped with a chair. Those systems work so that the user gives a final destination and supervises as the wheelchair is in charge of the complete navigation task [START_REF] Wang | Hybrid map-based navigation for intelligent wheelchair[END_REF]. In many cases, the user expresses explicitly his intended destination by pointing it into a map or choosing it from a menu. This type of interaction is not very natural and sometimes can be complicated for the user (specially elder users).

Semi-autonomous systems use the estimation of the user's plan as a key point because it allows the automatic controller to adjust its actions to the desire of its user. Some methods to perform an implicit estimation of the user's intention from joystick commands have been proposed in [START_REF] Demeester | Bayesian Estimation of Wheelchair Driver Intents: Modeling Intents as Geometric Paths Tracked by the Driver[END_REF], [START_REF] Demeester | Bayesian plan recognition for Brain-Computer Interfaces[END_REF]. They model the user's intention as a set of possible trajectories. A probability distribution is maintained over the set of trajectories and the most likely one is selected within a Bayesian framework.

In [START_REF] Taha | POMDP-based longterm user intention prediction for wheelchair navigation[END_REF] a learned Partially Observable Markov Decision Process (POMDP) was used to estimate the user's intended destination in a topological map of the environment. Places of interest are selected as those locations in the environment where the user spends comparatively most of his time. The user drives the wheelchair from one spatial location to another while the robotic device avoids obstacles in the middle.

Other smart wheelchairs limit their assistance to collision avoidance and leave the majority of planning and navigation duties to the user. These systems do not normally require prior knowledge of an area or any specific alterations to the environment. They require instead more planning and continuous effort on the part of the user and are only appropriate for users who can effectively plan and execute a path to a destination [START_REF] Vanacker | Context-based filtering for assisted brain-actuated wheelchair driving[END_REF].

B. Human Aware Navigation

Hall [START_REF] Hall | The hidden Dimension: Man's Use of Space in Public and Private[END_REF] classified human interactions based on a concept of distance, creating the concept of "public distance" to refer to situations in which people give a speech, and "social distance" that characterizes situations in which people talk to each other for the first time. Our approach uses the notion of social distances to find an appropriate location to place the wheelchair when the user wants to interact with the other people in a respectful and comfortable way.

[10] proposes a model of approaching behavior with which a robot can initiate a conversation with people who are walking. To prevent failures, their model includes a prediction of the walking behavior of people, choosing a target person, planning its approaching path, and nonverbally indicating its intention to initiate a conversation.

Other work, [START_REF] Hansen | Adaptive human aware navigation based on motion pattern analysis[END_REF], introduced an adaptive system which detects whether a person seeks to interact with the robot based on the person's pose and position, that system was presented as a basis for human aware navigation. Their results showed that the system was capable of navigating according to past interaction experiences and to adapt to different behaviors.

A method for detecting and tracking groups of people was presented in [START_REF] Luber | Multi-Hypothesis Social Grouping and Tracking for Mobile Robots[END_REF]. They address the problem of detecting and learning socio-spatial relations between individuals and Fig. 1. System Architecture Overview tracking their group formations using multiple social grouping hypotheses in a recursive way. This method also improves person-level tracking in two ways: the social grouping information is fed back to predict human motion over learned intra-group constraints and to support data association by adapting track-specific occlusion probabilities.

Personal space, o-space and their relation to comfort were addressed in [START_REF] Rios-Martinez | Understanding human interaction for probabilistic autonomous navigation using Risk-RRT approach[END_REF] where a risk based navigation was extended to include risk due to discomfort. Human movement is supposed to be known by learning of typical trajectories in a particular environment.

III. SYSTEM ARCHITECTURE

Fig. 1 presents an overview of the complete system. Several subsystems were developed by our team (shown in orange) while the other necessary modules were taken from the open source community (white blocks).

Destination Inference [START_REF] Escobedo Cabello | Multimodal Control of a Robotic Wheelchair: Using Contextual Information for Usability Improvement[END_REF]: This subsystem estimates the desired destination within the map of the environment among a list of possible predefined goals and meeting points computed by the social filter. The probability for each destination is computed using a Dynamic Bayesian Network that considers the current position of the wheelchair, user's face direction and the initial voice command. The destination with the highest probability is selected and sent to the navigation module.

Social Filter [START_REF] Rios-Martinez | Understanding human interaction for probabilistic autonomous navigation using Risk-RRT approach[END_REF]: Detects social interactions and creates virtual obstacles corresponding to those interaction zones in order to produce socially acceptable motion. The social filter integrates constraints inspired by social conventions to evaluate the risk of disturbance and takes it into account when performing the autonomous navigation planning.

Social Grid: Receives information from the social filter and projects it to an occupancy grid adequate to be used by standard motion planning algorithms.

Meeting Points: The social filter uses geometrical information from groups of people to compute good places for the wheelchair to be placed in order to join the group.

A. Destination Inference

The user's intention is modeled as a set of possible destinations. Those destinations may be defined by the user's habits (places where the user spends most of his time during the day), interesting points taken from the map of the environment as doors, desks and other facilities and dynamic destinations related to humans in the scene -For more details consider reading our previous work [START_REF] Escobedo Cabello | Multimodal Control of a Robotic Wheelchair: Using Contextual Information for Usability Improvement[END_REF]. In this work we propose a method to detect and use destinations defined by social groups of people.

The presented reasoning method is based on a dynamic Bayesian network described in Fig. 2.

The joint probability distribution that corresponds to the graph (Fig. 2) is:

P (D 0:t , C 0:t , X 0:t ) = P (X 0 ) P (V 0 ) P (D 0 |X 0 , V 0 ) * ∏ k=1:t [P (D k |D k-1 )P (C k |D k , X k ) P (X k )] (1)
Where, X t is the position of the wheelchair in a 2D plane with respect to the map reference frame. This position is discrete and its domain is defined by a grid whose size is according to the map of the environment. It is assumed that the position of the wheelchair X t is a fully observable measurement given by the localization system.

V 0 is the starting vocal command. This variable can take two possible values {GO, JOIN}. C t represents the direction pointed by the user at time t. This direction is centered on the user's local frame.

The variable D represents both the static destinations and the meeting points. The domain of the D variable (denoted as Ω D ) is defined as the union between the set of static destinations Ω S and the set of meeting points Ω M .

Ω D = Ω S ∨ Ω M (2)
Where, Ω D is the domain of the D variable. Ω s = {d (1) , d (2) , ..., d (N ) } is the set of N static destinations. Ω m = {d (N +1) , d (N +2) , ..., d (N +M ) } is the set of M meeting points.

The inference question to be answer at time t is P (D t = d (i) |C 1:t , X 0:t , V 0 ). It can be recursively computed Using Bayes rule as:

P (D (i) t |C 01t , X 0:t , V 0 ) = αP (X 0:t ) P (V 0 ) * P (D 0 |X 0 , V 0 ) P (C t |X t , D (i) t ) * ∑ j [P (D (i) t |D (j) t-1 )P (D (j) t-1 |C 1:t-1 , X 0:t-1 , V 0 )] (3)
The term P (D 0 |X 0 , V 0 ) will be referred as the initial probability model, the term

P (C t |X t , D (i) t ) is the command model, P (D (i) t |D (j) t-1
) is the transition model. The initial probability model has two different modalities, when V 0 = GO the value of this term is encoded in a probability table [START_REF] Escobedo Cabello | Multimodal Control of a Robotic Wheelchair: Using Contextual Information for Usability Improvement[END_REF], whose values are obtained by learning typical destinations in the environment according to the user's habits. However , when V 0 = JOIN the meeting points get an evenly distributed high probability while all the others get low values.

P (D (i) 0 |X 0 , V 0 = GO) = P (D (i) 0 |X 0 ) (4)
P (D (i) 0 |X 0 , V 0 = JOIN ) = { α/M if d (i) ∈ Ω m (1 -α)/N otherwise (5)
In our tests we used α = 0.9.

P (C t |X t , D (i)
t ) represents the probability that a command C t will be sent by the user when he is located at position X t and his destination is d (i) at current time t.

P (C t |X t , D (i) t ) = P ′ (C t |X t , D (i) t ) ∑ i P ′ (C t |X t , D (i) t ) (6)
where,

P ′ (C t |X t , D (i) t ) = 1 σ √ 2π exp -1 2 ( a i σ ) 2 (7)
The a i term is the angle between the command and the destination (the value of the standard deviation σ was set from experimental data to have a value of 0.2).

P (D

(i) t |D (j)
t-1 ) expresses the probability that the current estimated destination changes respect the last one.

P (D (i) t |D (j) t-1 ) = { (gain -1)/(N + gain -1) if D (i) t = D (j) t-1 1/(N + gain -1) otherwise (8)
Where, N is the number of possible goals in the environment. gain expresses how important the previous estimation will be with respect to the new data. Large values of this term will lead to slow responses whenever the user changes his intended destination. After practical experimentation it was defined a value of P (D

(i) t |D (j) t-1
) that is 10 times bigger if the last estimated goal was the same.

The selection of the goal is performed using exact inference. The selected destination is the one with highest probability.

d sel = arg max d (i) P (D t = d (i) |C t , X t ) (9)
This goal is then sent to the autonomous navigation system when the user gives one of the vocal commands "GO" or "JOIN".

B. Social conventions in human navigation

In normal situations, people try to maintain a minimum space between them to feel comfortable. Hall [START_REF] Hall | The hidden Dimension: Man's Use of Space in Public and Private[END_REF] modeled the level of discomfort experienced by the person related to the distance at which he is approached. This idea was formalized as the concept of personal space, which characterizes the space around a human being in terms of comfort during social activities. Similarly, interactions between two or more people have been modeled using the concept of o-space by [START_REF] Kendon | Spacing and orientation in co-present interaction[END_REF].

1) Personal Space: The implemented model consists in blending two Gaussian functions Γ f and Γ b (Eq.10) both of them centered in the position of the person. The first one represents the personal space in front of a human and is wider than the last one representing the back space. The Gaussian values represent the risk of disturbance associated with a point in the space around the pedestrian.

Γ x,Σ (p) = e -1 2 (p-x) t Σ -1 (p-x) (10)
where x, p are in R 2 and Σ is a diagonal covariance matrix defined as:

Σ = ( σ 2 x 0 0 σ 2 y ) . ( 11
)
The appropriate shape of the personal space is obtained by selecting the same values for σ x in both Γ f and Γ b but different values for σ y , being the one of Γ f the double of the value for Γ b .

2) O-Space: When two or more people are in a conversation, they tend to adopt formations with circular shapes. The o-space could be taken as a circle whose center coincides with that of the inner space. For the specific case of two people, some formations, called F-formations, have been identified as being particularly frequent [START_REF] Kendon | Spacing and orientation in co-present interaction[END_REF]. The social filter identifies individual F-formations (Vis-a-vis or V-Shape) and builds the corresponding o-space. In Fig. 3, the calculated ospace for a Vis-a-Vis interaction is shown. The o-space is represented by a two-dimensional Gaussian function Γ c of covariance matrix S and centered in C, then for each point Q around the center we have:

Γ C,S (Q) = e -1 2 (Q-C) t S -1 (Q-C) (12)
where S is a diagonal covariance matrix as defined in Eq.11. To get the shape of the O-space depending on the formations, some values has been chosen for the parameters which are shown in the next table:

Formation σ x σ y Vis-a-vis D H /3 2 * DH/3 V-Shape D H /3 2D i C. Social Cost Grid
Many planning algorithms use occupancy grids as its basic data structure to represent the state of the world. Practical implementations of such algorithms typically used in robot navigation systems (in this work A* and dynamic window) consider occupancy grids.

In order to adequate the data computed by the social filter we project the values of the personal space and interaction space into a 2D occupancy grid. Two different thresholds are used as shown in Fig. 4. Cells in the occupancy grid with personal and/or interaction space values higher than the threshold are marked as occupied zones so the planner will avoid passing through those regions, while regions with lower values are marked as free space.

Occupied regions computed in the social grid are then combined with the occupancy grid created with the robot's laser range finder and the map of the environment. In order to consider the size of the wheelchair the regions marked as occupied are inflated by an 'inflation radius' according to the size of the robotic wheelchair (Fig. 4(b)). Threshold applied to the Gaussian functions delivered by the social filter to created the occupancy-cost grid. (b) The inflated area (green) considers the dimensions of the wheelchair. In order to guarantee a path free of collisions the center of the wheelchair should never be inside this area. Fig. 5. Projection of both p-space and o-space on the occupancy grid used by the planning algorithm. The Gaussian in the middle represents the value of the interaction space, which is projected on the occupancy grid as occupied space (blue area) and then inflated to consider the size of the robot (green area).

D. Computing Meeting Points

A geometrical model is used to determine those places where the robotic wheelchair should be placed to attract the attention and become part of an interacting group. Those points are located in an area, roughly coincident with the ospace of an interaction where the robot can share the space in an equitable way with the other humans already present.

Fig. 6 shows different types of social formations and where the meeting points should be located. Fig. 6(a) represents a vis-vis or frontal formation where two humans H 1 and H 2 are talking face to face. In this case two meeting points are placed on the line perpendicular to the H 12 line that joins H 1 and H 2 . the points are located at a distance D H /2 from the center of the o-space. D H is the distance between the two persons.

In the case of a v-shape formation Fig. 6(b) a meeting point is placed on the line formed by V i and H 12 . V i is the intersection of the two lines of sight of interacting humans. The computed meeting points are used by the inference algorithm as possible destinations as it will be described later. IV. SIMULATION EXPERIMENTS The proposed approach was designed considering the experimental scenario shown in Fig. 8. People in the scene were tracked to learn the typical destinations that are then placed in the map of the environment (red circles in Fig. 8(a)). Each destination has a related probability value as described in section III-A. The values for the presented experimental scenario where presented in a previous publication [START_REF] Escobedo Cabello | Multimodal Control of a Robotic Wheelchair: Using Contextual Information for Usability Improvement[END_REF].

The user can start the movement at any location of the experimental scenario, he drives the wheelchair by seeing towards his desired destination and saying a vocal command as JOIN or GO.

A. Face Control Subsystem

The user can control the robotic wheelchair by using the movements of his face. The direction of the face is estimated by a random forest classifier which takes as input a 3D point cloud (for further information on this subject please refer to [START_REF] Fanelli | Real time head pose estimation with random regression forests[END_REF]).

B. Voice Control Subsystem

As it is impossible to deduce the full range of user intentions only analyzing the user's face direction; it is necessary the system to be multi-modal. For this reason, we used the Google voice recognition service with a reduced set of instructions -related to our application-to give those commands that it would be difficult to express using only the face as: stop, move, faster, etc. Some words have an strong context related to the type of destination that they can be applied to, for example the verb GO is more general than the verb JOIN which means that it does not give so much information to infer the user's desired destination. However, when the user uses the keyword JOIN, it is more likely that he intends to approach a group of people if there is any present on the scene. This behavior is encoded in the probability distribution of Eq. 5.

C. Multimodal Interface

In the example (Fig. 9), the user is looking in the direction of the arrow (blue) and there are two people in the middle of the simulated INRIA-hall. The task of the wheelchair is to navigate towards the goal that gets the highest posterior probability.

When the user gives the vocal command to start the movement (GO or JOIN) the prior probabilities are loaded and the posterior probability is computed considering the direction of his face, the probability for each goal is depicted in Fig. 9 as the size of the sphere. The goal with the highest posterior probability is sent to the navigation module, the wheelchair plans the path and compute the necessary 9. In this example we show the different modalities that can be achieved by our method. In the first example (a); given the strong context related to the use of the verb JOIN the probabilities are computed in such a way that those destinations defined by a group of people get a higher value. However the use of the verb GO, (b); does not give such contextual information; therefore, the destination that gets the maximum probability is mainly defined by the direction in which the user is pointing and the prior probability model. velocity commands to arrive to the destination while avoiding obstacles and respecting the personal and interaction spaces. The destination inference method considers both meeting points defined by the two persons in vis-vis formation. In the first case Fig. 9(a): the user gives a JOIN order. Given the strong context related to this verb the probabilities are computed in such a way that the meeting points get a higher probability value than all the other possible destinations.

In the second case Fig 9(b) the user gives a GO order which does not give as much contextual information as the JOIN command, so the probability distribution is more evenly distributed, therefore, the destination that gets the maximum probability value is mainly defined by the direction in which the user is looking.

In case of an error, for example when the inferred destination can not be reached because it is occupied or there is not enough evidence to decide which destination is better. The system will ask for help and wait for the user to give a new command. In the example Fig. 10 the system selects a meeting point that is not reachable by the wheelchair, then it asks for help and waits for the user to give a new command. When the new command is considered a new destination could be selected and if the new destination is valid the wheelchair will start moving.

V. CONCLUSIONS AND FUTURE WORK

The approach presented in this paper integrates a human aware cost grid with an standard autonomous navigation controller. The system was designed to improve both usability taking advantage of the user intention inference method and sociability by including the concepts of personal space and o-space in the navigation system.

The main contributions presented in this article were; first of all, the method to compute socially accepted points to join a group of people and then use them as possible destinations for the user. Second, the use of information from several user-machine interfaces -voice and face position-, to reduce ambiguities in the inference model.

Fig. 2 .

 2 Fig. 2. Dynamic Bayesian Network used to estimate the user intended destination Dt.

Fig. 3 .

 3 Fig. 3. O-space calculated by the Social Filter Module for a Vis-a-Vis formation. The maximum risk of disturbance is located at o-space center, in the picture the disturbance is represented by the Gaussian function height.

  Fig. 4.Threshold applied to the Gaussian functions delivered by the social filter to created the occupancy-cost grid. (b) The inflated area (green) considers the dimensions of the wheelchair. In order to guarantee a path free of collisions the center of the wheelchair should never be inside this area.

Fig. 6 .

 6 Fig. 6. Geometrical placement of meeting points in different social formations. (a) For a vis-vis formation. (b) For a v-shape formation. The location of Meeting points depends on the focus of attention defined by the orientation of the bodies φ.

Fig. 7 .

 7 Fig. 7. The resulting meeting points for a vis-vis formation (black markers). If the robot reaches that position, a group of three will be formed. (a) Lateral view. (b) Top view.

Fig. 8 .

 8 Fig. 8. Possible experimental scenario is located at INRIA main entry hall where typical destinations have been detected in previous works [2]. (a) The lines in the figure represent some of the trajectories normally followed by people in the scene. Typical destinations (end of a trajectory) are marked with circles. (b) In the intended application, the position of the user's face and the vocal command are used to infer his desired destination. Social formations are detected and the corresponding meeting points (d(10), d(11)) are considered as possible destinations by the destination inference algorithm.

  Fig.[START_REF] Hall | The hidden Dimension: Man's Use of Space in Public and Private[END_REF]. In this example we show the different modalities that can be achieved by our method. In the first example (a); given the strong context related to the use of the verb JOIN the probabilities are computed in such a way that those destinations defined by a group of people get a higher value. However the use of the verb GO, (b); does not give such contextual information; therefore, the destination that gets the maximum probability is mainly defined by the direction in which the user is pointing and the prior probability model.

Fig. 10 .

 10 Fig. 10. Error handling. When the system detects an error in the inferred destination (it is blocked or there is not enough evidence to decide) (a) it asks for help to the user (b). The user gives a new command (c) and if the new destination is valid it will move towards it.

This work has been funded by the INRIA large scale initiative project PAL, Personally Assisted Living