Chloé Mimeau

Georges-Henri Cottet

Iraj Mortazavi
email: iraj.mortazavi@cnam.fr

G.-H Cottet

Direct numerical simulations of three-dimensional flows past obstacles with a vortex penalization method

Keywords: flow past a sphere, flow past a hemisphere, vortex methods, particle methods, penalization method, semi-Lagrangian methods

Direct numerical simulations of three-dimensional flows

past obstacles with a vortex penalization method

Introduction

Vortex methods are Lagrangian methods consisting in discretizing the vorticity on particles which move with the flow and carry local vorticity values. They approximate the vorticity formulation of the Navier-Stokes equations (Vorticity Transport Equation) focusing the computational task on the vortical zones, thus allowing to reduce the simulation time. They are ideal for vortex-dominated flow simulations, like flows past bluff-bodies, airfoils, propellers, turbine profiles as well as complex internal configurations, thanks to a lower numerical diffusion and the dynamic adaptivity of the particles to the support of the vorticity. Moreover, continuity and far-field boundary conditions are satisfied by construction.

Vortex methods have been widely used to study wakes and more generally advection dominated flows and several works have provided these methods a solid mathematical framework to understand their convergence [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF].

For external flows, vortex methods often rely on the use of Fast Multipole Methods (FMM). This technique, introduced by [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], allows to reduce the algorithm complexity of the computation of the velocity of the vortex elements from O(N 2 ), obtained with a classical resolution of the Biot-Savart law, to O(N log(N )) . This has been widely used in the context of vortex methods [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF][START_REF] Marzouk | K-means clustering for optimal partitioning and dynamic load balancing of parallel hierarchical N-body simulations[END_REF][START_REF] Rossinelli | Multiresolution flow simulations on multi/many-core architectures[END_REF][START_REF] Yokota | FMM-based vortex method for simulation of isotropic turbulence on GPUs compared with a spectral method[END_REF]. Alternatively, vortex methods can rely on grid-based Poisson solver to compute velocity fields [START_REF] Ould-Salihi | Blending finite-difference and vortex methods for incompressible flow computations[END_REF]. In that case vorticity values have to be interpolated to the grid and velocity values have to be interpolated back to particles. In particular [START_REF] Cottet | Advances in direct numerical simulation of 3D wall-bounded flows by Vortex-In-Cell methods[END_REF] indicates comparisons of CPU times between FMM and methods based on Poisson solvers which indicate that, except if the support of vorticity has very small dimension compared to the computational box, Poisson-based vortex methods are more economical in 3D than FMM-based vortex methods, by several order of magnitudes. However [START_REF] Yokota | FMM-based vortex method for simulation of isotropic turbulence on GPUs compared with a spectral method[END_REF] shows that, when implemented on GPU with very large numbers of particles, FMM-based methods can close the gap with vortex methods based on Poisson solvers.

Vortex methods, and more generally particle methods, also face accuracy issues related to particle distortions. One popular method to handle this issue is to remesh frequently -in general at each time-step -particles on regular grids. This method was used in [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] to obtain reference results for flows past a 2D circular cylinder. With this technique, particle methods become semi-Lagrangian methods and their accuracy can be rigorously analyzed in terms of the interpolation kernels used to remesh particles [START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] under stability conditions on the time-step which do not rely on the interparticle spacing. In remeshed vortex methods, vorticity grid values obtained after particle remeshing can be readily used to obtain the velocity through a grid-based Poisson solver. Alternative methods allow to adapt velocity and diffusion kernels to the distribution of particles [START_REF] Schrader | Discretization correction of general integral PSE operators for particle methods[END_REF][START_REF] Rossi | The diffused vortex hydrodynamics method[END_REF] or to adapt the particle distribution to the boundaries [START_REF] Rossi | Simulating 2D viscous flow around geometries with vertices through the Diffused Vortex Hydrodynamics method[END_REF] (see also [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF] for earlier references on this issue) . However, to our knowledge, this type of methods has not been implemented yet in 3D flows.

A third issue faced by vortex methods, and more generally by all methods based on the vorticity formulation of the Navier-Stokes equations, is related to boundary conditions. The no-slip condition can be turned into a vorticity flux boundary condition that can be handled together with the vorticity diffusion by integral equations [START_REF] Koumoutsakos | Viscous boundary conditions for vortex methods[END_REF][START_REF] Graziani | From a boundary integral formulation to a vortex method for viscous flows[END_REF]. This method has been used with success in a number of two and three-dimensional simulations using vortex particles [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF][START_REF] Cottet | Advances in direct numerical simulation of 3D wall-bounded flows by Vortex-In-Cell methods[END_REF]. These simulations use body-fitted particle distributions. Alternative methods, primarily defined and used in the context of grid-based Eulerian methods [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Fadlun | Combined immersedboundary finite-difference methods for three-dimensional complex flow simulations[END_REF], consist in viewing obstacles as immersed boundaries which appear as a forcing term in the Navier-Stokes equations. Immersed boundary methods can easily be implemented on a non-conforming grid or particle distribution, albeit at the expense of a lower accuracy near the boundary. The Brinkman penalization method is a particular case of immersed boundary methods where the velocity is penalized inside the obstacle. Penalization methods can be interpreted as originally formulated in velocity-pressure formulation and can be extended to vorticity formulations [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkman penalization[END_REF] and to vortex methods [START_REF] Coquerelle | A vortex level-set method for the twoway coupling of an incompressible fluid with colliding rigid bodies[END_REF].

Vortex methods using a penalization method to enforce the no-slip condition have been developed and used in the context of fluid-structure interaction problems [START_REF] Coquerelle | A vortex level-set method for the twoway coupling of an incompressible fluid with colliding rigid bodies[END_REF][START_REF] Gazzola | Simulations of single and multiple swimmers with non-divergence free deforming geometries[END_REF]. In [START_REF] Rossinelli | GPU accelerated simulations of bluff body flows using vortex particle methods[END_REF] the algorithmic simplicity of the penalization method was utilized to enable an efficient GPU implementation of the method for 2D bluff-body flows. In [START_REF] Mimeau | Passive Flow Control Around a Semi-Circular Cylinder Using Porous Coatings[END_REF][START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF] the method was in particular used to model flows around turbines or to optimize porous layers for passive control to reduce drag behind a semi-circular cylinder.

In the present method we study flows past a sphere and an hemisphere and we specifically focus on two new aspects of the method which significantly reduce the computational cost of the method: we propose a simple method relying on FFT-based Poisson solvers to approximate the far field boundary conditions and we extend the directional remeshing technique introduced in [START_REF] Magni | Accurate, non-oscillatory remeshing schemes for particle methods[END_REF][START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] to the three-dimensional Navier-Stokes equation.

This paper is organized as follows. In the first three sections we recall the governing equations, we indicate how the body forces are computed in such a model and we detail our treatment of the far-field boundary conditions. In section 5 we detail the algorithm and the numerical schemes to discretize the problem. Section 6 is devoted to the numerical results. We perform a validation and a convergence study of the proposed method for threedimensional flows past a sphere. We then present similar investigations in the context of flows around a hemisphere. In this section we also discuss the computational cost and parallel scalability of the method. Finally, some conclusions and perspectives related to the present work are outlined in the last section.

Governing equations

This study is based on the vorticity formulation of the incompressible Navier-Stokes equations in a domain D, which read:

∂ω ∂t + (u • ∇)ω -(ω • ∇)u = 1 Re ∆ω in D. (1)
In this equation ω, u and Re respectively denote the vorticity, the velocity and the Reynolds number. This equation has to be coupled to the system giving the velocity in terms of the vorticity. Using the incompressibility condition, the velocity may be directly linked to the vorticity through the following Poisson equation:

∆u = -∇ × ω. (2)
The system (1)-( 2) has to be complemented by appropriate boundary conditions, both at solid boundaries and, if necessary, at artificial boundaries. This latter issue will be addressed in detail in section 4. The modeling of incompressible flow around an obstacle is realized in this work using the Brinkman penalization method [START_REF] Caltagirone | Sur l'intéraction fluide-milieu poreux : Application au calcul des efforts exercés sur un obstacle par un fluide visqueux[END_REF][START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. This technique relies on the simple following idea: the fluid in which the flow evolves is assumed to be a porous medium of infinite permeability while, on the other hand, the solid obstacles immersed in the fluid are considered as media with zero permeability. A flow evolving in such media may be modeled by the Brinkman-Navier-Stokes equations, which stands for the whole domain and which contains an additional term in the Navier-Stokes equations, acting as a forcing term. The latter, called the penalization term, directly comes from the Darcy equations that govern flows in porous media. In this work, the penalization term is expressed using vorticity formulation. One therefore obtains the non-dimensional penalized Vorticity-Transport-Equations (or Brinkman-Navier-Stokes equations), originally proposed by [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkman penalization[END_REF]:

∂ω ∂t + (u • ∇)ω -(ω • ∇)u = ∇ × λ χ b (u b -u) + 1 Re ∆ω in D, (3) ∆u
= -∇ × ω in D, (4)
where χ b denotes the characteristic function that yields 0 in the fluid and 1 in the solid body, u b indicates the rigid body velocity and λ = µΦH/ρku ∞ is the non-dimensional penalization parameter, with k the intrinsic permeability, µ the viscosity, Φ the porosity of the porous material, H the height of the obstacle, ρ the fluid density and u ∞ the main uniform fluid flow velocity.

The main advantage of this method is that it needs neither the mesh to fit the boundaries nor to specify no-slip boundary conditions. Moreover, the penalized Vorticity-Transport-Equation (3) appears to be very convenient to model the flow in the whole domain thanks to the dimensionless penalization factor λ, whose value allows to distinguish between the different materials. In this study we set H = ρ = u ∞ = 1 and we recall that the porosity Φ is close to 1 as imposed by the Brinkman equation [START_REF] Nield | Convection in Porous Media[END_REF]. Therefore λ essentially depends, in the inverse proportion, on the intrinsic permeability k of the medium. Varying the value of λ thus directly defines the different media. Indeed, in the fluid, the intrinsic permeability coefficient k goes to infinity, thus the fluid can be considered numerically as a porous media with a very high permeability. We set λ = 0 in this region. As a consequence, the penalization term vanishes in equation ( 3), and we naturally recover the dimensionless Vorticity Transport Equation (1). On the contrary, the solid has a permeability coefficient k close to zero, it can be consequently modeled by fixing the penalization parameter λ to a very high value. In this study λ is set to 10 8 in the solid. It was proved in [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] that solving equation [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] with such a value of λ was equivalent to solve Darcy's law in the solid. As a conclusion, at a given flow regime, the variation of λ corresponds to the variation of k and specifies the intrinsic porous material permeability. Figure 1 gives a schematic representation of the overall physical representation and the associated numerical settings. The physical description (left picture) involves the physical quantities H, u ∞ , Φ, k and ν = µ/ρ. The associated numerical setting (right picture) involves the characteristic function χ b enabling to geometrically define the solid body as well as the penalization term λ = νΦH/ku ∞ obtained from the adimensionalization of the penalized Navier-Stokes equations and allowing to numerically specify the permeability of the considered medium. For further details about the derivation of the penalized Vorticity-Transport-Equations from the Brinkman equation, the reader is referred to the work [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF].

Forces evaluation

The accurate computation of the unsteady aerodynamic force exerted by the fluid on an immersed body is important in the framework of this study. Indeed, the force provides a reliable diagnostic to validate the ability of the proposed method to model the problem correctly. The method adopted in this work to compute the aerodynamic force in three-dimensions is the numerical approach introduced by [START_REF] Morency | Aerodynamic force evaluation for ice shedding phenomenon using vortex in cell scheme, penalisation and level set approaches[END_REF] under the name of change of momentum. It is based on the fact that the time dependent force exerted by a fluid on a body B of constant mass is directly proportional to the body acceleration. For unit density one has:

F = - d dt B u dx. (5)
If we think in terms of numerical iterative process, we can consider that, at each time step, the fluid is allowed to enter inside the rigid body before the penalization term enforces the velocity u to be equal to u λ inside the body. Therefore, the obstacle experiences the following change of momentum ∆q:

∆q = B (u -u λ ) dx, (6)
and the force is finally obtained by the formula:

F = d dt B (u -u λ ) dx. (7)
The drag and lift coefficients are then defined by:

C D = F x Q × A , C L = F y Q × A , C S = F z Q × A , (8)
where the lift coefficient is decomposed into C L (vertical lift) and C S (side lift), Q = 1 2 ρu 2 ∞ denotes the dynamic pressure and A is the reference area. As this force calculation method is inherently based on the penalization technique, that is to say on the continuous velocity enforcement inside the body, it appears as a natural choice in the framework of the present numerical method. However, it may show some limits when high Reynolds numbers are considered, due to the grid refinements needed to accurately take into account the obstacle boundaries. In this case, other formulations may be used, like for instance the one proposed by [START_REF] Noca | A comparison of methods for evaluating time-dependent fluid dynamic forces on bodies, using only velocity fields and their derivatives[END_REF] where the forces are evaluated on a control volume that includes the solid body.

Far-field boundary conditions

In the present study we consider flows past obstacles in a domain with a prescribed inlet normal velocity and periodic boundary conditions in the transverse directions. To obtain optimal computational efficiency we will use periodic boundary conditions in all directions, which will allow to rely on simple FFT to obtain the velocity from the vorticity. This however requires appropriate corrections to account for the inflow and outflow boundary conditions.

Solenoidal eddies absorption

A first consequence of periodic boundary conditions is the emergence, in the upstream region, of eddies coming periodically from the outlet (we suppose in this study that the domain has a sufficient size in the y and z directions so that the wake do not reach the concerned boundaries). In order to discard these arising vortices, one needs to apply a specific treatment on the vorticity field. This treatment consists in an absorption, in the flow direction, of the vorticity within a band located at the outlet of the domain.

More precisely we proceed as follows. Given a periodic vorticity field ω periodic and its associated velocity field u periodic we first define:

ũ = f (x) u periodic + (1 -f (x)) u ∞ , ( 9
)
where f is the following one-dimensional smoothing function defined in the flow direction [START_REF] Parker | Stability and vortex shedding of bluff body arrays[END_REF]:

f (x) =          1 if x < x b tanh(α(x -x c )) -tanh(α(x e -x c )) tanh(α(x b -x c )) -tanh(α(x e -x c )) if x b ≤ x ≤ x e 0 if x > x e . (10)
In the above expression x b , x c and x e respectively refer to the beginning, the center and the end of the absorption band at the outlet (see Figure 2). The parameter α allows to adjust the steepness of the absorption function f . One then obtains a corrected vorticity field:

ω = ∇ × ũ (11) = f (x) ω periodic + [∇f (x) × u periodic ] + ∇[1 -f (x)] × u ∞ (12) =   ω x f (x) ω y f (x) ω z f (x)   +   0 -f (x)u z f (x)u y   +   0 f (x)u z∞ -f (x)u y∞   . ( 13
)
Note that the above absorption guarantees the solenoidal condition div ω = 0, a property which would not be satisfied if the absorption mechanism [START_REF] Cottet | Advances in direct numerical simulation of 3D wall-bounded flows by Vortex-In-Cell methods[END_REF] was applied directly to the vorticity.

Velocity correction

Our choice of periodic boundary conditions in the Poisson solver to compute the velocity from the vorticity requires to adjust the correct flux through the computational box. Moreover, periodic boundary conditions also imply a zero circulation, i.e. D ω = 0. This condition is not satisfied as soon as eddies leave the computational box and a correction is needed to account for the non-zero circulation. In the sequel x denotes the direction along which the inlet velocity is prescribed and y and z the transverse directions where periodic boundary conditions are assumed.

To derive the correction enforced on the velocity calculations we introduce the following notation:

u = u + u. (14)
In this relation, u denotes the output velocity field obtained from the FFTevaluations performed in the Fast Poisson Solver and u represents the spaceaverage velocity, that we aim to correct. The x-axis of the Cartesian system corresponds to flow direction and the velocity is initialized setting

u(t = 0) = (u x∞ , u y∞ , u z∞ ).
Since u x is a spatially constant flux, then ∂u x ∂x = ∂u x ∂y = ∂u x ∂z = 0 and the mean vorticity is given by:

ω x = ∂u z ∂y - ∂u y ∂z , ω y = - ∂u z ∂x , ω z = ∂u y ∂x . ( 15
)
First of all, a correction of the streamwise component of the velocity u x is necessary in order to impose the correct flux at the inlet. This correction is determined as follows.

Correction for the streamwise velocity u x

Let I be the inlet rectangular surface of the parallelepipedic computational box, normal to the flow direction, of area L y L z . L y and L z respectively denote the lengths of the computational box in the y and z directions. According to the decomposition u = u + u and knowing that u x is a constant and that the desired flow rate at the inlet is equal to u x∞ L y L z , one has:

I u x dydz desired flow rate = I u x dydz calculated flow rate + I u x dydz, (16)
so

u x = u x∞ -I u x dydz L y L z , (17)
and finally:

u x = u x∞ + u x -I u x dydz L y L z . ( 18
)
Next we need to correct the transverse velocities u y and u z in order to recover a non-zero circulation.

Correction for the spanwise velocity u y

According to relations (15) one has u y = ω z x + c 1 , with c 1 a constant value that we aim to evaluate and

ω z = 1 L x L y L z D ω z dx, with L x L y L z
corresponding to the volume of the computational domain D. From the decomposition u = u + u, one therefore concludes that u y = u y + ω z x + c 1 . Thus:

I u y dydz = I u y dydz + I ω z x dydz + I c 1 dydz. ( 19
)
Since the flow rate is desired to be u y∞ L y L z , then:

c 1 = u y∞ -I u y dydz L y L z -ω z x 0 , (20)
and finally the corrected y-component of the velocity fields stands as:

u y = u y∞ + u y + ω z (x -x 0 ) -I u y dydz L y L z , (21)
where x 0 corresponds to the x-coordinate of the inlet surface I in the Cartesian system.

Correction for the transverse velocity u z

Following the same principle, and according to the relations (15) one obtains u z = -ω y x + c 2 with c 2 a constant value. From the decomposition u = u + u, one gets u z = u z -ω y x + c 2 . Finally, mimicking the process used for u y , one gets:

u z = u z∞ + u z -ω y (x -x 0 ) -I u z dydz L y L z . ( 22
)
In order to validate the use of periodic boundary conditions combined with the above corrections, we consider the flow past a 3D sphere at Re = 300 in a bounded domain D = [-2, 8.24] × [-2.56, 2.56] × [-2.56, 2.56]. The target free stream velocity is uniform at the inlet, with u ∞ = (1, 0, 0). Figure 3 shows upstream isocontours of u x at the end of the simulation (T = 75) in the inlet YZ plane (x = x 0 = -2) for two different cases. In Figure 3 we have plotted the values of the streamwise velocity obtained at the inlet. In the first case (left picture), the absorption band lies from x b = 7.24 to x e = 8.24 and α = 10 in the definition of f (eq. 10). In this case the values of u x lie between 0.973 and 1.009. In the second case (right picture), the absorption band is twice thicker (x b = 6.24, x e = 8.24) and u x inlet ∈ [0.977, 1.007]. The actual inlet streamwise velocity thus only shows a maximum relative error of respectively 2.7% and 2.3% for the first and the second case compared to the target uniform free stream flow (u x∞ = 1). Let us also notice the symmetrical patterns of the velocity isocontours. This symmetry is explained by the periodic boundary conditions. Simulations were also performed with a domain twice longer in the x-direction, i.e. x b = 17.48, x e = 18.48 and very similar results were obtained (u x inlet ∈ [0.975, 1.004]). Finally, a simulation was carried out in the original domain D with a higher Reynolds number, Re = 1000, setting x b = 7.24, x e = 8.24 and comparable results were achieved (u x inlet ∈ [0.977, 1.004]). From this study one can conclude that, in view of the speed up obtained by this method, the proposed treatment of the inlet and outlet boundary conditions give satisfactory results.

The vortex penalization method

Numerical discretization

In this section we present the vortex penalization method used to discretize the 3D penalized Vorticity-Transport-Equations ( 3)-( 4). In this study, we will only consider the case of fixed obstacles, which means u b = 0. The method can be extended to the case of moving bodies along the lines of [START_REF] Coquerelle | A vortex level-set method for the twoway coupling of an incompressible fluid with colliding rigid bodies[END_REF][START_REF] Gazzola | Simulations of single and multiple swimmers with non-divergence free deforming geometries[END_REF]. Equations ( 3)-( 4) are discretized on a parallelepipedic computational box. The domain is meshed with a uniform Cartesian grid and periodic boundary conditions are prescribed on the walls as described in section 4. These equations are discretized through a viscous splitting algorithm, consisting in solving successively the different operators within the same time iteration [START_REF] Chorin | Numerical study of slighly viscous flow[END_REF]. Concerning the main steps of this splitting algorithm, the advection is performed in a semi-Lagrangian way through a vortex method with particle remeshing at each time step to avoid the distortion effects [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. The penalization equation, the stretching and the diffusion terms as well as the Poisson equation are solved on the underlying grid. Algorithm 1 gives the successive discretization steps.

In this work, the Poisson equation (step [START_REF] Morency | Aerodynamic force evaluation for ice shedding phenomenon using vortex in cell scheme, penalisation and level set approaches[END_REF] in algorithm 1) is solved on the grid using Fast Fourier Transforms (FFT) and periodic boundary conditions. Following section 4, the discretization of this equation is preceded by a vorticity absorption (step (26)) and followed by a velocity correction (step [START_REF] Noca | A comparison of methods for evaluating time-dependent fluid dynamic forces on bodies, using only velocity fields and their derivatives[END_REF]). The diffusion equation (step [START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF]) is also discretized using FFTbased evaluations. In the following, we will focus on the specific features that contribute to the efficiency of this method for 3D bluff body flows.

Algorithm 1 Initialization

Preliminary velocity field (on grid)

u 0 = u ∞ = (1, 0, 0) (23)
Preliminary vorticity field (on grid) ω 0 = (0, 0, 0)

Vorticity penalization (on grid)

ω 0 λ = ω 0 + ∇ × -λχ∆tu 0 1 + λχ∆t (25)
Marching Vorticity absorption (on grid)

ω n = ∇ × [f u n + (1 -f )u ∞ ] (26) Poisson problem (on grid) ∆u n = -∇ × ω n (27) Velocity correction (on grid) u n+1 = u n + u n (28)
Force computation (on grid)

F = d dt B (u n+1 -u n+1 λ ) dx (29)
Vorticity penalization (on grid)

ω n λ = ω n + ∇ × -λχ∆tu n+1 1 + λχ∆t (30)
Stretching (on grid)

∂ω n ∂t = div(ω n : u n+1 ) (31)
Diffusion (on grid)

∂ω n ∂t = 1 Re ∆ω n (32)
Advection (on particles)

∂ω n+1 ∂t + (u n+1 • ∇)ω n+1 = 0 (33)
Evaluation of ∆t adapt (on grid) ∆t

n+1 adapt = LCFL ∇u n+1 ∞ ( 34
)
Force computation (step ( 29))

As explained in section 3, the aerodynamic force is evaluated in 3D using the change of momentum equation [START_REF] Yokota | FMM-based vortex method for simulation of isotropic turbulence on GPUs compared with a spectral method[END_REF]. For this purpose, equation ( 5) is discretized using a 1 st Euler scheme:

F = - 1 ∆t B u n+1 -u n h 3 , (35)
and u n+1 is forced to be equal to the velocity inside the body through an implicit penalization scheme (cf equation [START_REF] Kim | Laminar flow past a hemisphere[END_REF] in next subsection):

F = 1 ∆t B u n - u n 1 + λ∆t h 3 = B λu n 1 + λ∆t h 3 , (36)
where u n corresponds to the non-penalized velocity and h denotes the uniform grid space. This force expression is evaluated on the grid with a midpoint rule. We note that the above equation is consistent with the formulation [START_REF] Yokota | FMM-based vortex method for simulation of isotropic turbulence on GPUs compared with a spectral method[END_REF]: since in our calculations we have λ = 10 8 and ∆t ≥ 10 -3 , then λ∆t 1 and (36) gives:

F = B u n ∆t h 3 . ( 37
)
This is a 1 st order discretization of the change of momentum equation [START_REF] Yokota | FMM-based vortex method for simulation of isotropic turbulence on GPUs compared with a spectral method[END_REF] where u λ = 0. Let us also notice that, if ∆t → 0, then one would have:

F = B λu n h 3 , (38)
which corresponds to the discretized formulation of the force expression proposed in [START_REF] Caltagirone | Sur l'intéraction fluide-milieu poreux : Application au calcul des efforts exercés sur un obstacle par un fluide visqueux[END_REF] and obtained by integrating the penalization term on the volume of the body:

F = 1 Re B ∆u dx - B ∇p dx ≈ D λχ b u dx. (39)
In the present work one chooses equation [START_REF] Kim | Laminar flow past a sphere rotating in the streamwise direction[END_REF]. The force evaluation is performed immediately before the penalization step in algorithm 1 (step (29)).

Penalization (step (30))

In its velocity-vorticity formulation, the continuous penalization equation reads :

∂ω ∂t = -∇ × (λ χ b u). (40)
In order to discretize equation [START_REF] Campregher | Computations of the flow past a still sphere at moderate Reynolds numbers using an immersed boundary method[END_REF], we first discretize the penalization equation for the velocity, ∂ t u = -λχ b u. Following [START_REF] Coquerelle | A vortex level-set method for the twoway coupling of an incompressible fluid with colliding rigid bodies[END_REF] this is done using an implicit 1 st order Euler scheme:

u n+1 = u n 1 + λχ b ∆t . ( 41
)
This yieds the penalized velocity field that we denote by u λ . We can rewrite the above equation as:

u n+1 = u n - λχ b ∆t u n 1 + λχ b ∆t , (42)
which gives the following conservative expression for the penalized vorticity:

ω n+1 = ω n -∇ × λχ b ∆t u n 1 + λχ b ∆t . (43)
A 4 th order centered finite-differences scheme is used for the discretization of the curl operator. We notice that the velocity field is not directly penalized in this algorithm. The velocity penalization is implicitly realized through the resolution of the Poisson problem in step ( 27) of algorithm 1.

Stretching (step (31))

Several formulations may be used to compute the stretching term [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF]. In this study we choose the following conservative formulation:

∂ω ∂t = div(ω : u), (44)
where the notation ω : u stands for the tensor of component ω j u i . The time integration scheme chosen here to discretize this equation is the commonly used 3 rd order Runge-Kutta TVD (Total Variation Diminishing, i.e. non extra oscillations) scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. With this time discretization, the velocity field involved in the divergence operator is not modified. The divergence operator is discretized through a 4 th order centered finite-differences scheme on the grid.

Particle advection and remeshing procedure (step ( 33))

The semi-Lagrangian advection in the 3D case is performed using a 1 st order directional splitting approach [START_REF] Magni | Accurate, non-oscillatory remeshing schemes for particle methods[END_REF][START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF]. This technique consists in solving the convection/remeshing problem direction by direction. In other words one has to solve here three different mono-dimensional problems.

First, one solves the following differential equations to update the new positions of particles at time t n+1 from the grid positions x n i and the first component u 1 of velocity field u n at time

t n :    dx p dt = u 1 (x p ), t ∈ [t n , t n+1 ] x n p = x n . (45)
A 2 nd order Runge-Kutta scheme is used for the resolution of this equation. Once all the particles have been pushed in the same direction, they are remeshed on the grid, in this direction. Then one proceeds to the advectionremeshing in the second and third direction. The kernel chosen here is Λ 4,2 [START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF][START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF]. Its 1D support includes 6 points, it satisfies 4 moment conditions and is twice differentiable. It yields a second order approximation to the advection equation. It is defined by the following formulas: where |x| denotes the distance between the particle and a given grid point.

Λ 4,2 (x) =                    1 - 5
From this formula, one can evaluate direction-by-direction the weights collected by each grid point involved in the kernel support. In each direction, the remeshing scheme yields :

ω n+1 i = p ω n p Λ 4,2 x n+1 p -x i h . ( 46
)
This particle convection/remeshing algorithm allows a significant reduction of the computational efforts compared to the more traditional remeshing based on tensor product formulas.

Adaptive time step (step (34))

In particle methods, the non-linear stability of the convection scheme is ensured under the following condition on the time step:

∆t adv ≤ LCFL ∇u ∞ , (47)
where the Lagrangian CFL must satisfy LCFL ≤ 1. As rigorously analyzed in the convergence proof [START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF], this Lagrangian CFL is the appropriate condition under which consistency and stability of the remeshed particle method can be guaranteed. This Lagrangian stability condition imposes that particles trajectories do not cross. One can notice here that the time step in not constrained by the grid size or the distance between the particles but only by the flow strain, which often provides larger time steps compared to Eulerian grid-based schemes. In the present study, the time-step is updated at each iteration using the condition (47) with LCFL = 1/8. Except at the initial stage of the simulations, when strong vorticity layers develop around the obstacle, we observed that the time-step rapidly settles at a roughly constant value.

Implementation and library description

The vortex penalization method is implemented in a library which uses object oriented programming techniques in order to reach a high level of modularity, with a strong focus on usability and flexibility. Our goal is to enable the user to launch indifferently sequential or parallel (MPI) simulations. We use Python as an abstraction framework and the high level of abstraction provided by this language allows to conceal from the user the parallel paradigms and the low level implementations of the numerical algorithms. Concretely, the user only needs to describe his problem with high level components, mathematical operators, problem variables and physical domain description, and then has to choose the available numerical methods and algorithms, or develop his own modules, according to his needs. This abstraction framework enables us to develop a library that is portable to various kinds of modern computer architectures. More information about this library can be found in [START_REF] Cottet | High order semilagrangian particles for transport equations: numerical analysis and implementation issues[END_REF], where a special attention is made on the capability of the presented software to perform hybrid CPU-GPU numerical simulations. In the following, the computations are only CPU-based.

Numerical results for bluff body flows

In this section we present the numerical results obtained with the present method. Its convergence and accuracy are numerically analyzed for threedimensional bluff body flows. The numerical results presented in this section are based on DNS simulations. Laminar and transitional flows are considered for two distinct obstacles, namely the sphere and the hemisphere. In the hemisphere case, the bluff-surface faces the upstream flow and the flat one is oriented downstream. In both cases the non-dimensional diameter d of the obstacle and the free stream velocity are set equal to 1. The obstacles are centered at the origin of the computational box.

Numerical setup

In all our simulations the size of the computational box was chosen large enough so that the presence of the artificial boundaries was found to be negligible. A careful convergence study with respect to the domain size has been performed in [START_REF] Cottet | Advances in direct numerical simulation of 3D wall-bounded flows by Vortex-In-Cell methods[END_REF] for flow past a 2D circular cylinder using vortex methods. It demonstrates that a domain extending to about 6 cylinder diameters downstream is sufficient to provide reliable flow diagnostics. In the present work, the additional zone used for the absorbing boundary condition accounts only for about 10% of the total size (and thus of the total computational cost). In the subsequent simulations, the size of the domain is set to D = [-2, 8.24] × [-2.56, 2.56] × [-2.56, 2.56]. The computational domain D and the geometrical setup are shown in Figure 4 for both obstacles. The whole domain is discretized with a uniform Cartesian grid. For each of the different flows considered in this section, a convergence study will be carried out and comparisons of our results with respect to numerical and experimental references will be given.

Flow past a 3D sphere

Flow past a sphere is a common benchmark to validate the accuracy of a numerical method and to prove its capability to correctly model 3D bluff body flows.

Grid convergence study

A grid convergence study is performed in the context of flow past a 3D sphere at Re = 300. It involves four mesh sizes: h = 0.08, h = 0.04, h = 0.02 and h = 0.01, which correspond to a number of grid points ranging from 128 × 64 × 64 to 1024 × 512 × 512 in the domain D. Convergence orders are determined by computing the L 1 , L 2 and L ∞ norms of the error of two different flow quantities, namely the enstrophy Z = D |ω| 2 dx and the drag coefficient C D , with respect to the best resolved case (h = 0.01):

e(t) = || q best resolved (t) -q(t) ||, (48)
where the quantity q either denotes the drag coefficient C D or the enstrophy Z. In the particular context of this grid-convergence study, the time-step had to decrease proportionally to the grid size. Instead of (47), we used the following criterion to evaluate it:

∆t n adapt = CFL • h max 1≤i≤3 |u i | , with CFL = 0.5. ( 49
)
Figure 5 indicates that the method exhibits around second order convergence regarding drag coefficient (1.95 for L 1 , 2.46 for L 2 and 1.98 for L ∞ ) and between first and second order convergence for the enstrophy (0.93 for L 1 , 1.35 for L 2 and 1.06 for L ∞ ).

The asymptotic mean values of the drag coefficient and the enstrophy are given by Table 1 for each mesh size. Based on Table 1, one can consider that the grid convergence is roughly achieved for h ≤ 0.04. The next simulations of incompressible flows past a 3D sphere at Re = 300 will be performed with a mesh size set to 0.02 in order to better capture the boundary layer.

Re

Flow analysis and validations

Flow past a 3D sphere has been widely studied numerically and experimentally which allowed to precisely determine the different flow regimes. Following the studies of [START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF] and the description given in [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF], we recall the regimes related to Reynolds numbers ranging from 0 to a few thousands: • Re < 210 -212: the flow is steady and axisymmetric. An axisymmetric separation bubble exists at the back of the sphere with zero lift force.

• 210 -212 < Re < 270 -290: the flow is steady and planar-symmetric.

A pair of planar-symmetric vortical structures appears in the wake.

• 270 -290 < Re < 350 -375: the flow is unsteady, time periodic and planar-symmetric. The wake is characterized by a cyclic shedding of ring vortices, with planar-symmetry.

• 350 -375 < Re < 800: the flow is unsteady, non-periodic and fully asymmetric.

• Re > 800: the flow is unsteady, non-periodic, asymmetric and a Kelvin-Helmholtz instability occurs in the shear layer, propagating in the wake. The latter becomes fully turbulent, characterized by the emergence of small scales.

In the sequel, we propose validation studies for two Reynolds numbers corresponding to clearly distinct regimes: Re = 300 and Re = 1000.

Flow past a sphere at Re = 300

We first present the results obtained at Re = 300. The simulation setup is the following: the 3D computational box D = [-2, 8.24] × [-2.56, 2.56] × [-2.56, 2.56] is meshed by an uniform 512 × 256 × 256 Cartesian grid (corresponding to h = 0.02). As in [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF], a perturbation is addressed to the uniform upstream flow between the non-dimensional time T = 3 and T = 4, in order to trigger the instability. This perturbation is applied on the y component of the velocity and stands as u y∞ = sin(π(T -3)). The time step is defined according to equation (47). It is approximately equal to 0.011 all along the simulation. The penalization parameter λ is set to 10 8 .

Table 2 compares the values of the drag and vertical lift coefficients as well as the Strouhal number obtained using the present method with results in literature. [START_REF] Roos | Some experimental results on sphere and disk drag[END_REF] 0.629 --Johnson & Patel [START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF] 0.656 -0.069 0.137 Tomboulides & Orszag [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] 0.671 -0.136 Constantinescu & Squires [START_REF] Constantinescu | LES and DES investigations of turbulent flow over a sphere[END_REF] 0.655 -0.065 0.136 Kim & Choi [START_REF] Kim | Laminar flow past a sphere rotating in the streamwise direction[END_REF] 0.657 -0.067 0.134 Ploumhans et al. [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF] 0.683 -0.061 0.135 Present work 0.673 -0.066 0.133 As the table shows, our results coincide well with these references. Figure 6 depicts the ω z and ω x isocontours for every quarter period of a shedding cycle. The contours of the streamwise (ω x ) and spanwise (ω z ) components of the vorticity obtained with the present method are compared with those published by [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF]. As can be observed on the figure, the vortical features are in very good agreement with the reference results. Moreover, one can notice that a true periodic regime is attained (the snapshots in Figures 6a and6b depict the vorticity field at five equispaced times within one shedding period) and that the wake is perfectly symmetric in the XZ plane (see the ω x isocontours in the XZ plane in Figure 6b). These observations confirm the description of the flow regime corresponding to Reynolds 300 thoroughly investigated by One can notice in particular that, at Re = 300, the vortex shedding from a sphere does not symmetrically alternate along the flow direction (XY plane), where the top eddies are much stronger than the bottom ones (Figure 7a). This observation confirms the one made by [START_REF] Kim | Laminar flow past a sphere rotating in the streamwise direction[END_REF] and corresponds to the fact that the vertical lift coefficient CL is non-zero (see Table 2). Flow past a sphere at Re = 1000

Re = 300 Authors cD cL S t Roos & Willmarth
The flow analysis at Re = 1000 is now considered. The simulation is performed on an uniform 1024 × 512 × 512 Cartesian grid (corresponding to h = 0.01). A perturbation is applied between T = 3 and T = 4 on the y component of the velocity, defined by u y∞ = 0.1 sin(π(T -3)). The adaptive time step returns a value roughly equal to ∆t = 0.005, which corresponds to a CFL number equal to 0.6. The penalization parameter λ is fixed to 10 8 . A plot of the time average streamwise velocity u x along the centerline is given in Figure 8a. The results obtained with the present method are compared to those of [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] together with the experimental data of [START_REF] Wu | Sphere wakes in still surroundings at intermediate Reynolds numbers[END_REF] at Re = 960. It can be seen that the length of the recirculation zone, characterized by negative u x values, obtained in the present case (1.5 diameters) is similar to the one obtained in [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] (approximately 1.7 diameters). Furthermore, the u x values for x/d ≥ 2 are in a very good agreement with other numerical and experimental studies ( [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] and [START_REF] Wu | Sphere wakes in still surroundings at intermediate Reynolds numbers[END_REF]) and show the accuracy of the present method. The three experimental and numerical results are however slightly different for smaller x/d values. As the number of points in the boundary layer are almost the same for both numerical methods (3.3 points in the present case and 3.5 in [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF]) these discrepencies may be due to the interpolation errors in the boundary layer that contains very high velocity gradients. Moreover, some experimental errors at x/d < 2 also could be induced during the measurements.

Concerning the time evolution of force coefficients, Figure 8b shows that the time evolution of the drag coefficient C D and the vertical lift coefficient C L coincide with the one found by [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF]. As emphasized in [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF], although the prescribed perturbation is symmetric, the flow quickly looses its plane symmetry. This is confirmed by Figure 9, giving the evolution of the force coefficients for a larger time range. The side lift force coefficient C S is nonzero and shows significant variations which are actually in the same order of magnitude as C L . Concerning the drag coefficient C D , the mean value obtained in the present study is 0.485, which agrees well with the numerical values reported by [START_REF] Poon | Numerical studies of high Reynolds number flow past a stationary and rotating sphere[END_REF] and [START_REF] Campregher | Computations of the flow past a still sphere at moderate Reynolds numbers using an immersed boundary method[END_REF], respectively equal to 0.46 and 0.478. Isocontours of |ω| in the near wake are depicted in the XY and XZ plane for early times, respectively in Figure 10a and10b. They are compared to those of [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF]. These figures, which show a good agreement, reveal that the flow is initially symmetric in the XZ plane. However, as highlighted by Figure 11, the near wake becomes asymmetric at T = 18. The same value of T is reported by [START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF]. This loss in planar-symmetry increases in time and is perfectly clear at T = 40 (see Figures 11 and12). These observations on the wake symmetry are directly related to the C S time evolution given in Figure 9 where, from

Flow past a 3D hemisphere

Hemispherical bodies correspond to configurations that may be used in several engineering applications. In a previous work [START_REF] Mimeau | Passive Flow Control Around a Semi-Circular Cylinder Using Porous Coatings[END_REF], this geometry was used to investigate drag reduction around the side view mirror of a ground vehicle. Due to the presence of a flat back wall with sharp edges, the flow past a hemisphere is a steeper problem compared to the flow past a sphere and is therefore challenging from a numerical point of view. Vortex methods adapted to geometry singularities have been proposed in two dimensions in [START_REF] Rossi | Simulating 2D viscous flow around geometries with vertices through the Diffused Vortex Hydrodynamics method[END_REF]. The purpose of the present work is in particular to measure to which extent penalization methods can handle this type of singularity in three dimensions, at low and moderate Reynolds numbers. The refinement study proposed in the following section proves that the accuracy of the penalization method for sharp eddy geometries is the same as the one obtained for smooth bluff bodies.

Grid convergence study

A grid convergence study is performed for flow past a 3D hemisphere at Re = 300. As for the sphere, four different mesh sizes are considered for this convergence study and the time step is defined according to equation (49). The order of convergence for the drag coefficient is found to be 1.80 for L 1 , 2.23 for L 2 , 1.48 for L ∞ (Figure 13a) and 0.95 for L 1 , 1.49 for L 2 , 1.23 for L ∞ concerning the enstrophy (Figure 13b). These convergence orders are very similar to those found for the sphere in the previous section. They are complemented by Table 3, giving the converged mean values of the two flow quantities at Re = 300. Here again, on the basis of these results one may consider that the grid convergence is roughly obtained by setting h = 0.04.

Flow analysis and validations

To our knowledge there are only few experimental and numerical results devoted to the problem of flow past a 3D hemisphere. This section is based on the study carried out by Kim & Choi [START_REF] Kim | Laminar flow past a hemisphere[END_REF], who investigated in details the behavior of flows past a hemisphere at different Reynolds numbers, ranging from 100 to 300. In their study, the authors report the following flow regimes:

• Re = 100, 130, 150, 170: the flow is steady and axisymmetric.

• Re = 180, 190: the flow is steady and planar-symmetric.

• Re = 200, 210: the flow is unsteady, time periodic and asymmetric.

• Re = 220, 250, 270: the flow is unsteady, time periodic and planarsymmetric.

• Re = 280, 300: the flow is unsteady, non-periodic and asymmetric.

Let us first draw some comparisons between these flow regimes and the ones of the sphere, listed in the previous section. One can notice that steady flows over a hemisphere are very similar to those over a sphere, even if the transitions between the different flow regimes occur at lower Reynolds numbers in the hemisphere case. Indeed, as explained in [START_REF] Kim | Laminar flow past a hemisphere[END_REF], the steady axisymmetric flow is characterized by an axisymmetric recirculation zone behind the hemisphere, while for the steady planar-symmetric flow two streamwise vortical structures are formed in the wake. However, according to the authors, an important difference exists between the unsteady flows past a hemisphere and those past a sphere. We indeed notice that when 220 ≤ Re ≤ 270, the flow past a hemisphere recovers planar-symmetry while for lower Reynolds number (Re = 200, 210), it was totally asymmetric. This behavior is not observed in the case of the sphere.

Flow past a hemisphere at low Reynolds numbers

In this validation study we will focus on three values of the Reynolds number, corresponding to different unsteady flow behaviors: Re = 200, 250, and 300. These Reynolds numbers are also studied in details in [START_REF] Kim | Laminar flow past a hemisphere[END_REF]. For each simulation, we consider a 512×256×256 Cartesian grid (corresponding to h = 0.02). As for the sphere, the flow instability is triggered by a perturbation imposed between T = 3 and T = 4.

The flow at Re = 200 is first considered. Figure 14a shows the time evolution of the drag and lift coefficients, which are compared to those obtained by [START_REF] Kim | Laminar flow past a hemisphere[END_REF]. First, one can notice a very good agreement between the different results, which is confirmed by the first lines of Table 4 giving the mean values of the force coefficients at Re = 200. Concerning the Strouhal number, one can see in our case and in the reference results that the Strouhal based on the lift is half the one based on the drag. Finally, we focus on the side lift coefficient C S . At Re = 200 the time average of this coefficient is not zero, which corroborates the fact that the flow does not maintain the planar-symmetry. This statement is confirmed by Figure 14b the plane XZ at this Reynolds number.

Let us now focus on flow past a hemisphere at Re = 250. Figure 15a and Table 4 show that the present outcomes coincide well with the reference results. As for the flow at Re = 200, the shedding frequency is equal to the frequency of the lift (St lift = 0.128) and half the frequency of the drag (St drag = 0.256). Finally, Figure 15b confirms that the flow at Re = 250 recovers the symmetry in the XZ plane. This observation explains the zero time average side lift coefficient ( CS = 0) read in Table 4.

At Re = 300, according to Figure 16a the flow has a more complex behavior, showing drag and lift coefficients with complicated variations in time. The vorticity norm reported in Figure 16b also highlights the absence of any planar-symmetry, as expected. On the other hand, Table 4 reveals that, at this regime, there co-exist two distinct values of the Strouhal number based on the power spectrum of the drag. The lower frequency (St drag1 = 0.04) is due to the large-scale instability of the wake, while the higher frequency (St drag2 = 0.270) is attributed to the small-scale instabilities caused by the separation of the shear layer. The St drag2 value is still approximately equal to twice the dominant frequency of the lift (St lift = 0.138). One can also remark in Table 4 that the time average lift coefficient CL is zero (like for Re = 200 and Re = 250), which is not the case for flow past a sphere at Re = 300 (see CL in Table 2). This difference can be explained by the discrepancy in the process of vortex shedding. Indeed, as can be seen in Figure 16b, in the case of a hemisphere the vortices shed alternatively with the same strength in the flow direction (XY plane), leading to CL = 0. On the contrary, at the same Reynolds number, we shown in the previous section that the vortex shedding from a sphere was not alternating symmetrically (see Figure 7a).

Flow past a hemisphere at Re = 1000

Our last investigation is performed at Re = 1000. To to our knowledge, no study has been carried out experimentally or numerically at such regime, and more generally at a Reynolds number larger than Re = 300. On the basis of the grid refinement study presented in section 6.3.1 we have chosen for this case a grid-step value h = 0.01. All the other numerical parameters used for this simulation are the same as the one chosen for the sphere at the same Reynolds number (see section 6.2.2). The results presented in Figure 17 indicate that the characteristics of the flow past a hemisphere at Re = 1000 are similar to the one observed in the case of the sphere at the same regime, in the sense that it is unsteady, non-periodic and asymmetric. Figure 17a shows the streamwise velocity u x along the centerline, averaged on the time range T = [60, 100]. Although this mean velocity profile is evaluated on a later time range, one can nevertheless notice that it is comparable to the one of the sphere (see Figure 8a), with in particular a recirculation zone extending from the rear back of the obstacle to x 2. The time history of the force coefficients is reported in Figure 17b. It shows that from T 55 the wake becomes chaotic and is characterized by important and non-periodic variations of the side lift coefficient C S . Three-dimensional representations of the turbulent wake are given in Figure 17c at T = 60 and T = 100. Figure 17c also reports the time averages of the vorticity magnitude between T = 70 and T = 110 in the XY and XZ plane, highlighting the chaotic feature of the wake. In particular, we remark in both planar section on the right pictures of Figure 17c that a large recirculation area with low vorticity values exists behind the hemisphere, followed by an important zone characterized by high vorticity values.

We note here that, under the same initial conditions and identical numerical parameters, the wake of the hemisphere looses its symmetry in the XZ plane much later than the one of the sphere (T 18, see Figure 11). Figure 18 depicts the isocontours of the streamwise vorticity past a hemisphere in the XZ plane at different times. One can verify in particular that, as suggested by the evolution of the side lift coefficient C S (see Figure 17b), the wake looses its symmetry between T = 40 and T = 60. The retarded asymmetric behavior can be explained by the flat back wall of the hemisphere and seems triggered by the insight of secondary instabilities in the shedding shear areas.

Computational resources

The simulations presented in this paper were performed on a cluster made of bi-Xeon Sandy-Bridge cores. Table 5 gives the numerical setup, the resources and the CPU time-to-solution for different bluff body flows analyzed in this work. The simulations presented in this paper were realized in the domain D = [2, 8.24]×[-2.56, 2.56]×[-2.56, 2.56]. The different t end values reported in the table and corresponding to the non-dimensional final time of the simulations were chosen according to the reference test-cases we wanted to compare with. ∆t adapt denotes the mean of the non-dimensional adaptative time step. This mean value is calculated on the time range [ t, t end ], where t denotes the time from which the flow regime is established. For simulations performed at Re = 300, the allocated memory was less than 2.7 Gbytes and for simulations carried out at Re = 1000, it did not exceed 21 Gbytes. The CPU time-to-solution shown in this table correspond to a CPU time of about 40 seconds per time step on 64 cores for the high resolution cases with N = 1024 × 512 2 , and of about 8 seconds per time step on 32 cores for the lower resolution cases with N = 512 × 256 2 . This indicates a weak scalability in this range of resolutions around 80%. Figure 19 gives the percentage of the different stages implied in the resolution of the 3D penalized Vorticity-Transport-Equations ( 3)-(4), in one simulation time step for the execution on 1 core. In particular, one can notice that the FFT grid-based Poisson solver only represents 5% of the total CPU cost of one time step.

Conclusion

In this work, a vortex penalization technique has been employed to simulate three-dimensional bluff body flows. The method combines the simplicity of the penalization technique and the robustness of remeshed particle methods. The penalization approach allows the use of simple Cartesian grids and associated fast Poisson solvers to compute velocity fields. For this particular case of flows where the boundary conditions are either periodic or prescribed through upstream flows and free downstream flow, we showed that it is pos-T=20 T=40 T=60 T=80 T=100 at Re = 1000 was finally performed. The numerical results reported in this work prove the ability of the present method to correctly account for the expected behavior of the flow at transitional or highly transitional regimes and to accurately handle geometries showing sharp edges. A study of the computational cost of the method and of its scalability on parallel computers confirms that the method provides a versatile tool to analyze bluff-body flows.

However, when increasing the Reynolds number it can easily be foreseen that using uniform Cartesian grid would quickly require prohibitive mesh sizes in order to correctly capture the thin boundary layer in DNS. Domain decomposition techniques and/or adaptive particle methods, in the spirit of [START_REF] Bergdorf | Multilevel adaptive particle methods for convection-diffusion equations[END_REF][START_REF] Bergdorf | A Lagrangian particle-wavelet method[END_REF][START_REF] Rossinelli | Multi-resolution adapted grids for remeshed vortex methods on multicore architectures[END_REF], are appealing extensions of remeshed particle methods which will become mandatory in these regimes. The challenge will then be to reach a high degree of adaptivity while keeping the simplicity and associated good scalability of the penalized remeshed particle methods. Another, complementarity direction that we plan to investigate is the implementation in the context of vortex methods of second order penalization methods [START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF].

Figure 1 :

 1 Figure 1: Schematic representation of the physical description (left) and its numerical parametrization (right) for flow past a solid body using the penalization method framework.

Figure 2 :

 2 Figure 2: Schematic representation of the smooth vorticity absorption performed at the outlet.

Figure 3 :

 3 Figure 3: Flow past a 3D sphere: isocontours of the inlet streamwise velocity u x in the YZ plane when imposing periodic boundary conditions coupled with adapted corrections. (Left picture) The width of the absorption band equals 1. (Right picture) The width of the absorption band equals 2.

Figure 4 :

 4 Figure 4: Computational domain D.

Figure 5 :

 5 Figure 5: Grid-convergence study for flow past a sphere at Re = 300. The errors e L 1 (black), e L 2 (red) and e L ∞ (blue) are plotted against the grid step. (a) Convergence of the drag coefficient (b) Convergence of the enstrophy.

  (a) Isocontours of ω z in the XY plane (b) Isocontours of ω x in the XZ plane

Figure 6 :

 6 Figure 6: Vorticity isocontours for flow past a sphere at Re = 300 for every quarter period of a shedding cycle. Comparison between [4] (left columns) and present results (right columns).

[ 33 ,

 33 [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF]] and given at the beginning of this section. Finally, Figure7depicts 3D views of the vorticity norm |ω| resulting from simulations performed in a longer computational domain, D = [-2, 18.48]×[-2.56, 2.56]×[-2.56, 2.56].

Figure 7 :

 7 Figure 7: 3D instantaneous vorticity norm |ω| at T = 72 for flow past a sphere at Re = 300.

Figure 8 :Figure 9 :

 89 Figure8: Flow past a sphere at Re = 1000: (a) Average streamwise velocity u x along the x-axis: comparison of the present results (red curve) with numerical results[START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF] (green curve) and experimental data[START_REF] Wu | Sphere wakes in still surroundings at intermediate Reynolds numbers[END_REF] at Re = 960 (black circles). (b) C D (solid lines) and C L (dashed lines) time evolution: comparison of the present results (red curves) with[START_REF] Ploumhans | Vortex methods for direct numerical simulation of three-dimensional bluff body flows: Applications to the sphere at Re = 300, 500 and 1000[END_REF] (black curves).

  (a) in the XY plane (b) in the XZ plane

Figure 10 :Figure 11 :

 1011 Figure 10: Isocontours of |ω| for flow past a sphere at Re = 1000. Comparison between [4] (left columns) and present results (right columns) at T = 6, 8, 10, 12.

Figure 12 :

 12 Figure 12: 3D instantaneous vorticity norm |ω| at T = 40 for flow past a sphere at Re = 1000.

Figure 13 :

 13 Figure 13: Grid-convergence study for flow past a hemisphere at Re = 300. The errors e L 1 (black), e L 2 (red) and e L ∞ (blue) are plotted against the grid step. (a) Convergence of the drag coefficient. (b) Convergence of the enstrophy.

  Comparison of the time evolution of the force coefficients with the numerical results of[START_REF] Kim | Laminar flow past a hemisphere[END_REF]. Instantaneous 3D vortical structures with isocontour |ω| = 0.8.

Figure 14 :

 14 Figure 14: Flow past a hemisphere at Re = 200.

  Comparison of the time evolution of the force coefficients with the numerical results of[START_REF] Kim | Laminar flow past a hemisphere[END_REF]. Instantaneous 3D vortical structures with isocontour |ω| = 0.8.

Figure 15 :

 15 Figure 15: Flow past a hemisphere at Re = 250.

  Comparison of the time evolution of the force coefficients with the numerical results of[START_REF] Kim | Laminar flow past a hemisphere[END_REF]. Instantaneous 3D vortical structures with isocontour |ω| = 0.8.

Figure 16 :

 16 Figure 16: Flow past a hemisphere at Re = 300.

Figure 17 :

 17 Figure 17: Flow past a hemisphere at Re = 1000. (a) Average streamwise velocity u x along the x-axis. (b) Time evolution of the force coefficients. (c) Instantaneous vortical structures at T = 60 and T = 100 and time average of vorticity magnitude between T = 70 and T = 110 in the XY and XZ plane.

Figure 18 :

 18 Figure 18: Isocontours of ω x in the XZ plane for flow past a hemisphere at Re = 1000.

Table 1 :

 1 Convergence study for flow past a sphere at Re = 300. CD and Z respectively denote the mean values of drag coefficient and enstrophy.

			= 300
	Grid	cD	Z
	h = 0.08	0.732 64.8
	h = 0.04	0.679 66.5
	h = 0.02	0.673 67.7
	h = 0.01	0.676 68.9

Table 2 :

 2 Comparison of drag and lift coefficients and Strouhal number for flow past a sphere at Re = 300. Star notation ( ) refers to experimental results.

Table 3 :

 3 Convergence study for flow past a hemisphere at Re = 300. CD and Z respectively denote the mean values of drag coefficient and enstrophy.

Table 4 :

 4 depicting the vorticity norm |ω| at two different times, where one can see that the wake is not symmetric in Comparison of mean drag and lift coefficients as well as Strouhal numbers for flow past a hemisphere at Re = 200, 250 and 300.

	Re = 200	cD	cL	cS	St lift St drag
	Kim & Choi [41] 0.790	0.0 ± 0.024	0.049	0.128 0.256
	Present	0.805 0.003 ± 0.025	0.043	0.129 0.257
	Re = 250	cD	cL	cS	St lift St drag
	Kim & Choi [41] 0.742	0.0 ± 0.05	0.0	0.127 0.254
	Present	0.757	0.0 ± 0.05	0.0	0.128 0.256
	Re = 300	cD	cL	cS	St lift St drag1 St drag2
	Kim & Choi [41] 0.715	0.0 ± 0.062	0.0	0.135	0.04	0.270
	Present	0.729 -0.002 ± 0.063 -0.0016 0.134	0.04	0.270

sible to rely on simple and efficient FFT based Poisson solvers for the evaluation of the velocity. Another original feature presented in this work concerns the remeshing process. The particle redistribution is indeed performed direction by direction, which allows significant computational savings in 3D compared to classical tensorial approaches.

The validation studies reported in the second part of this paper were first carried out in the case of flow past a 3D sphere. Grid-refinement studies enabled us to identify a global 2 nd order of convergence for the method. Careful comparisons with several numerical and experimental reference works in the literature showed satisfactory agreements. Similar studies were then performed considering a 3D hemisphere. Validations were supplied for Reynolds numbers between 200 and 300. The simulation of flow past a hemisphere


Éric Auteur

Leclerc

Smart City and E-Governance in India. Mapping a New Urban Landscape

Keywords: Smart city, E-governance, Web 2.0, Networks, India MOTS CLÉS Ville intelligente, e-gouvernance, web 2.0, réseaux, Inde

The creation of a new greenfield capital in Andhra Pradesh, Amaravati, is an opportunity for its Chief Minister to test new foms of e-governance through the mobilisation of ICT tools. According to Naidu, this project of "People's Capital", Praja Rajadhani in Telugu, will enable them to build a more inclusive city. To achieve its goals, he mobilises all the resources of the web to make the project more transparent and promises to rely on citizen participation. In this article we will evaluate by a structural analysis of the Internet networks built by this government, if these affirmations have a chance to be realised during the construction of the capital. The challenge is to move from a top-down approach, from government to citizens, set up by the e-government procedures of the 1990s, to a bottom-up approach, from the citizen to his government promised by the web 2.0.

INTRODUCTION

Since the creation of Telangana, the 29th Indian state, the residual State of Andhra Pradesh (AP), has had a decade to create a new capital. It is a new greenfield city, Amaravati, that N. Chandrababu Naidu, the State Chief Minister, has chosen to build. Based on the global model of the smart and sustainable city, Naidu wants to build an inclusive city which it named Praja Rajadhani, the "People's Capital" in Telugu. To achieve this goal, he opted for a set of ICT tools to improve the transparency of the building process and ensure an efficient monitoring system. Taking into consideration Naidu's past achievements in planning and urban governance through the example of Cyberabad, it seems legitimate to question the political mantra "People's Capital" and the concept of smart city mobilised for the future capital of Amaravati. In this paper, we will discuss the concept of "digital city" along its influence on the equity and inclusion definition of an Indian metropolis in progress. Our objective is to visualise this new e-governance machinery in order to understand if ICT tools deliver their promises of efficient and participative democracy.

SMART CITY AND E-GOVERNANCE: THE ANDHRA PRADESH CASE

If we take the definition of Dameri and Cocchia (2013), "a smart city is a well-defined geographical area, in which high technologies such as ICT, logistic, energy production, and so on, cooperate to create benefits for citizens in terms of well-being, inclusion and participation, environmental quality, intelligent development; it is governed by a well-defined pool of subjects, able to state the rules and policy for the city government and development". In a shorter way, we could say that a smart city improves spatial justice (inclusion and participation) through ICT. In their definition, they conceive the smart city with a bottom-up approach, from citizen to government. This participatory conception of the city is a direct consequence of the web 2.0 and a transition from e-government (paperless) to e-governance where citizens are engaged in the management of the city.

This the is road taken by the government of the new Andhra Pradesh. Chandrababu Naidu has embarked on a digital agenda to create the new capital, Amaravati. He tried to harness the web 2.0 with online engagement with citizens through multiple digital and social media platforms like Facebook (facebook.com/prajarajadhani), a twitter account (with hastags #manaamaravati, #manarajadhani) and a channel on Youtube (youtube.com/prajarajadhani). Of course, the questions are: are all these new online tools delivering their promises? Will Amaravati be a more inclusive city?

THE E-GOVERNANCE NETWORK OF AMARAVATI

As Amaravati is for now a paper dream, to answer these research questions we will analyse the various digital initiatives set up by the Andhra Pradesh State in order to realise its goal of e-governance. To conduct this survey we analysed the network of web entities concern by the project of Amaravati and performed a structural study of this network to compare the respective influence of official websites and other web actors [START_REF] Diminescu | E-diasporas Atlas : exploration et cartographie des diasporas sur les réseaux numériques Paris[END_REF][START_REF] Leclerc | Topological and Content Analysis of the Cyberspace of the Indian Diaspora[END_REF]. Our objective is to understand the range of the e-governance policy of the AP government, to whom these websites are related?

Mobilising a "curation-oriented crawler": Hyphe

Inspired by the methodology of digital humanities we decided to rely on the web as a field of study [START_REF] Ackland | Web social sciences: Conceps, data and tools for social scientists in a digital age[END_REF][START_REF] Rogers | Digital methods[END_REF]. The main method to analyse a network of online entities is to crawl the web in order to discover all the hyperlinks between these entities. For this purpose, we used Hyphe, a free web-tool1 developed by the Sciences Po Paris Medialab [START_REF] Jacomy | Hyphe, a curation-oriented approach to webcrawling for the social sciences[END_REF] which has been tailored for social scientists. Its developers define it as a "curation-oriented crawler" because "it uses the concept of 'corpus' as a pivot to articulate web-mining constraints with social science methodology". In this case, we decided to explore all the web entities related to the construction of Amaravati starting from the official websites of the AP government.

We proceeded in three phases: firstly we built the corpus combining automatic crawl and a manual selection process of the websites; which resulted in a final corpus of 100 web entities and 222 links. Secondly, we spatialised this network of web entities in order to be able to interpret the graph. For this second phase, we used Gephi which is an open graph visualisation platform 2 that mobilises the Force-Atlas algorithm developed by Jacomy to draw our graph. Thirdly, we interpreted of the graph to understand whether the e-governance network built buy the Andhra Pradesh government around Amaravati is open or closed. The main interest of Hyphe is that it is not a crawler dedicated to one kind of web entities (website, blog, or social networks, Facebook or Twitter) but it is open to the diversity of the web.

Lessons from the corpus building process

The construction of the corpus goes through a selection process which is already interesting. We have decided to keep only web entities directly related to the new capital Amaravati. Which explains why we eliminated many websites on Buddhism, including temples in various places of the world named Amaravati. They are not directly related to the new capital, but to the small city Amaravathi which was the capital of the Satavahanas from the 3rd century BC to 3rd AD, a small town 15 kms NW of the new capital. We also put aside the numerous websites concerning archeological museums or art galleries related to the remains of the stupa of Amaravathi that are nowadays stored in the British Museum. From all websites dedicated to architecture and architects, we only kept those involved in the project of Amaravati: Japanese architect Maki, and British architect Rodgers.

During the selection process, Hyphe enables you to make a first analysis by creating categories that describe web entities. At this stage, we defined three variables to characterise the corpus. First variable: the type, with two categories, website or social network (in other words the web 1.0 and the web 2.0). The second variable: the producer of the web entities, with 4 categories: government, private companies, associations or individuals. The third variable: the geographic origin of the web entities, with 4 categories: Andhra Pradesh, Maharashtra, India, World. The town of Amravati in Maharashtra was kept in the corpus as an element of comparison with the capital of Andhra Pradesh.

A TOP-DOWN APPROACH OF E-GOVERNANCE

Two generations of web-entities

As described in the previous section, the spatialisation of the network provides a graph (fig. 1) that we can analyse. The size of the nodes that represent web entities are proportional to the number of in-coming links, and also the label of the web entity:ig circles therefore representweb entities with maximum authority in the network, while other entities refer to them. The second element that we can analyse on this graph are the clusters of web entities that constitute aggregates [START_REF] Ghitalla | La géographie des agrégats de documents sur le web[END_REF]. They are built by the algorithm Force-Atlas because they have more interconnection within the cluster than with the rest of the network.

On this first representation we identified 5 clusters. The main cluster at the centre of the graph is made of websites (blue color) among which we can identify two authorities: the government of Andhra Pradesh (ap.gov.in) and the central government (india.gov.in). On the periphery, we find 4 clusters of social network entities. The bigger is created by the Facebook account of the capital Amaravati, Praja Rajadhani. There are also two smaller clusters of Facebook accounts one for N. Chandrababu Naidu (Facebook AP CM) on the left, and one for "Amaravati Voice" (a media dedicated to the capital and its development) on the right. The last cluster is made of 2 gephi.org twitter accounts. Therefore each social network has its own cluster because they form groups of interlinked accounts. And there are also for Facebook separate clusters of account. Chief Minister's and "Amaravati Voice"'s look like hub because one parent account has out-going links with other accounts. All Facebook clusters show the main account with its variants in the different countries. Web 1.0 and web 2.0 are not so represented with, respectively, 29 and 70 entities, even though web 1.0 seems more diverse that social network clusters that aggregate the same kind of accounts.

A State-sponsored network

The second representation (fig. 2) takes into account the producer of web entities. 50% entities are produced by government agencies both for the central cluster of the web 1.0. But these are also active on the web 2.0 with accounts on Facebook for the capital and for the Chief minister. The second main producers are private companies (24%), situated on the periphery together with the "Amaravati voice" cluster, and various isolated web entities websites of the architect, and of Surbana Jurong, a Singapore company that produced the Master plan of the capital). These private entities are not connected with the main network. The third kind of producers are individuals (16%), only present through twitter accounts, mainly a sample of people re-twitting official news about the capital. Finally the fourth producers are associations (7%), who are either not very active on the web, or without links to official web entities (and were therefore not within our scope).

A local network

The third representation is based on the geographic variable, where from are the producers creating these web entities? Four places are located on this graph (fig. 3): Andhra Pradesh, Maharastra, India and the world. This network is massively from Andhra Pradesh (81%) and especially thanks to social networks. This completes the previous analysis underlining the diversity of the web 1.0, and the homogeneity of social networks. For the geographic origin, the central cluster has two levels of government web entities, Andhra Pradesh et India. We also find 3 websites for the Maharastrian town of Amravati, isolated, on the right of the graph. But there is also the world scale, with links with the World bank or the British agency of development (DFID). The analysis of this network demonstrates that it was initiated by the Andhra Pradesh government both on the web 1.0 through e-governance portal and on social networks with official accounts on Facebook and Twitter. But this network is not well connected with the web because private companies keep themselves outside, and the bulk of the links are with other public entities (Indian government or international institutions).

CONCLUSION

Evaluating the inclusiveness of a policy is a rather difficult task, and even more when you consider a not-yet-realised project. The construction of Amaravati officially started on the 30th of September 2017; it will not be possible to verify the promises of the Chief Minister before at least a decade. This is why we decided to focus on the ICT tools used by this government to foster e-governance. Beyond the current statistics about the number of hits on a website or the number of followers for a social network, the structural analysis of the network of web-entities provide more insightful conclusions. In our case, we can conclude that the e-governance system of the new government of Andhra Pradesh is a small word in a sociological sense. The websites or accounts on social networks are few and they are heavily related with one another. It demonstrates that although the government e-governance strategy has been very active on the web, the reach is very poor. It used the web 2.0 but in a top-down approach and considers these tools as new ways to advertise its action, which basically remains very centralised. They might have a more participative e-governance process with some of the applications dedicated to the citizens, like "Mana Amaravati" 3 ; it is too early to assess it, but it will be interesting to analyse the data collected.

Figure 1 .

 1 Figure 1. Amaravati's network: types of web entities (built with Hyphe)

Figure 3 .

 3 Figure 3. Amaravati's network: geographic origins (built with Hyphe)

hyphe.medialab.sciences-po.fr

THE AUTHOR

Éric Leclerc

Université de Lille -TVES eric.leclerc@univ-lille.fr