Chloé Mimeau

Federico Gallizio

Georges-Henri Cottet

Iraj Mortazavi
email: iraj.mortazavi@cnam.fr

G.-H Cottet

G.-H Cottet

Vortex penalization method for bluff body flows

Keywords: vortex methods, penalization technique, moving obstacles, porous medium

come

INTRODUCTION

Vorticity is an important characteristic of flows in many applications. It is a signature of wakes in bluff body flows and, as such, is related to aerodynamic performance of moving or fixed obstacles [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF]. The vorticity formulation of the Navier-Stokes equations is therefore a natural framework to address questions related to flow control. In these equations, for moderate and high Reynolds numbers the vorticity advection is predominant, which motivates the use of Lagrangian or semi Lagrangian schemes to discretize these equations. Particle methods belong to this category. They to the case of a rotating turbine blade and section 5 is devoted to the analysis of the flow around a semi-circular cylinder with a porous coating, where the characteristics of the flows are analyzed from the point of view of flow control. Finally section 6 summarizes our results and draws some conclusions.

NUMERICAL ALGORITHM

Vorticity formulation and remeshed particle method

The dynamics of an incompressible flow without external force fields is governed by the incompressible Navier-Stokes equations :

∇ • u = 0 in D ( 1
)
∂u ∂t + (u • ∇)u = -∇p + 1 Re ∆u in D, (2)
where D is the computational domain and Re is the Reynolds number.

Taking the curl of equation ( 2) and using equation [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF], we get the Helmholtz or the Vorticity Transport Equation (VTE) for the vorticity ω = ∇ × u :

∂ω ∂t + (u • ∇)ω -(ω • ∇)u = 1 Re ∆ω in D, (3)
This equation has to be coupled to the system giving the velocity in terms of the vorticity. It is convenient to use the Helmholtz decomposition :

u = ∇ × ψ + ∇ϕ. (4)
The stream function ψ and potential ϕ satisfy the following systems in D :

-∆ψ = ω, ∇ • ψ = 0, ∆ϕ = 0, [START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF] complemented with appropriate boundary conditions. The VTE does not explicitly contain the pressure term but the pressure can be recovered from the velocity field by integrating the following Poisson equation :

∆p = -∇ • (u • ∇)u, in D, (6)
which is derived as the divergence of equation [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF].

The previous equations are approximated using a remeshed vortex particle method [START_REF] Chorin | Numerical study of slighly viscous flow[END_REF], [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF]. In this approach, the vorticity is carried by particles. Particles follow the trajectories associated to the velocity field. Simultaneously the vorticity carried by the particles is updated to account for vortex stretching and diffusion. Particle remeshing is used to overcome accuracy loss due to local distortion of the flow map [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. In that case, particles are remeshed through interpolation on a regular mesh at each time step. This grid is used to compute the velocity through a Poisson solver. Details of the algorithm are given in the sequel.

In incompressible flows, the unique sources of the vorticity are boundaries. It is therefore important to capture in an accurate way vorticity boundary conditions. In classical implementations of vortex methods, the no-slip boundary condition is satisfied through the creation of vortex elements [START_REF] Chorin | Numerical study of slighly viscous flow[END_REF][START_REF] Chorin | Vortex sheet approximation of boundary layers[END_REF][START_REF] Creusé | Vortex simulation of active control strategies for transitional backwardfacing step flows[END_REF] or by updating particle strength to account for vorticity fluxes at the boundary [START_REF] Koumoutsakos | Boundary conditions for viscous vortex methods[END_REF][START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | Advances in direct numerical simulation of 3D wall-bounded flows by vortex-in-cell methods[END_REF]. The no-through flow boundary conditions is implemented together with the Poisson equation to determine stream functions and potential in [START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF].

These methods involve several parameters to enforce the no-slip conditions at the boundaries : a more straightforward alternative to these techniques can be derived from the penalization method.

Vortex penalization method

Penalization in velocity formulation

In this section we show how the penalization method can be used successfully to model the flow of an incompressible fluid around an obstacle. In the penalization technique the system is considered as a single flow, subject to the Navier-Stokes equation with a penalization term that enforces continuity at the solid-fluid interface and rigid motion inside the solid. The main advantage of this method is that it needs neither the mesh to fit the boundaries nor to specify no-slip boundary conditions. In addition it allows to compute the pressure as a continuous field on the whole domain including the solids, and the lift and drag coefficients by integrating the penalization term inside the solid bodies [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].

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

In a porous medium the velocity u is given by the Darcy law :

u = - k µΦ ∇p, (7)
where p is the pressure, k is the intrinsic permeability, µ is the dynamic viscosity of the fluid and Φ is the porosity of the medium. On the other hand, assuming that the Boussinesq hypothesis is satisfied for the fluid saturating the porous medium, we get the Brinkman's equation that includes the diffusion inside the solid and is only valid when the porosity coefficient is close to one. Here we set this coefficient equal to one and take the hypothesis of the viscosity equivalence for both Darcy and diffusive terms ( [START_REF] Bruneau | Passive control of the flow around a square cylinder using porous media[END_REF][START_REF] Nield | Convection in Porous Media[END_REF]) :

∇p = - µ k u + µ∆u. (8)
Complementing the Brinkman's equation with the inertial terms and the Dupuit-Forchheimer relation [START_REF] Nield | Convection in Porous Media[END_REF], we finally obtain the modified Forchheimer-Navier-Stokes equation ( [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Bruneau | Numerical modelling and passive flow control using porous media[END_REF]) :

ρ ∂ t u + ρ (u • ∇) u + ∇p = - µ k u + µ∆u (9)
where ρ is the density of the fluid.

Then, using the following usual substitutions including the mean velocity of the fluid u and the size of the body H :

u = u u ; x = x H ; t = t H u ; p = p ρ|u| 2 , (10)
we get the penalized non-dimensional Navier-Stokes equations or Brinkman-Navier-Stokes equations for the unknowns (u , p ) denoted (u, p) for the sake of simplicity :

∂ t u + (u • ∇)u - 1 Re ∆u + u K + ∇p = 0 in D T = D × (0, T ) (11)
∇ • u = 0 in D T . ( 12
)
the penalization coefficient that defines the permeability of each region. In the fluid the permeability coefficient goes to infinity, the penalization term vanishes and we solve the non-dimensional Navier-Stokes equations. In the solid body the permeability coefficient goes to zero and it has been shown in [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] that solving these equations corresponds to solve Darcy's law in the solid and that the velocity is proportional to K. One can thus write K = 1 λ χ S where χ S denotes the characteristic function of S and λ is a penalization parameter satisfying λ 1. For numerical simulations we typically set λ ≥ 10 8 . If the full domain including the solid body is defined as D = F ∪ S, one finally obtains the following penalized Navier-Stokes equations :

∂ t u + (u • ∇)u - 1 Re ∆u + λχ S u + ∇p = 0 in D ( 13
)
∇ • u = 0 in D. ( 14
)
In addition as the pressure is computed inside the solid body, it is shown 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] that one can compute the drag and lift forces by integrating the penalization term on the volume of the body :

F = 1 Re S ∆u dx - S ∇p dx ≈ D λ χ S u dx, (15)

Penalization in vorticity formulation and implementation with a vortex method

We first show how to extend the penalization method to the vorticity formulation of the Navier-

Stokes equations.

By differentiation of the velocity based penalization equation [START_REF] Van Rees | Optimal shapes for anguilliform swimmers at intermediate reynolds numbers[END_REF], we obtain the following non dimensional [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with brinkman penalization[END_REF][START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] vorticity-penalization equation :

∂ω ∂t + (u • ∇)ω -(ω • ∇)u = 1 Re ∆ω + ∇ × [λ χ S (u s -u)]. (16)
In the above equation we have generalized the boundary value problem to allow for a rigid body motion of the obstacle S with velocity u s , where u s is the rigid body velocity and χ S is the characteristic function of the solid. This system has to be complemented by the usual system giving the velocity in terms of the vorticity :

∇ • u = 0 in D; ∇×u = ω in D (17)
Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls and appropriate boundary conditions on the walls Γ D .

We describe here how the vorticity equation ( 16) is discretized in a vortex method. For simplicity we assume that the computational domain, including the solid body, is covered by a uniform Cartesian grid with mesh-size h.

The equation ( 16) is split in two sub-steps. At each time step, one successively solves the following equations :

∂ω ∂t = ∇× (λ χ S (u s -u)) (18)
∂ω ∂t + (u • ∇) ω -(ω • ∇)u = 1 Re ∆ω. (19)
We set t n = n∆t, where ∆t is the time-step. If u n , ω n denote velocity and vorticity values on the grid at time t n , following [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] the resolution of equation ( 18) is performed using an Euler implicit scheme, which yields :

ω n+1 = ∇× u n + λ∆tχ S u n s 1 + λ∆tχ S . (20)
The right hand side above is evaluated by centered finite differences on the grid.

To solve the advection-stretching-diffusion equation [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], we proceed as follows. Particles are initialized on grid points where vorticity value is above a given threshold (typically 10 -5 ). If x n p , ω n p denote the locations (grid points) and vorticity of particles (ω n p = ω n+1 (x p )), the equation ( 19) is firstly handled by solving the following system of ODEs [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF] :

       dx p /dt = u p dω p /dt = [∇u p ]ω p . ( 21
)
This system is typically solved using a second order Runge-Kutta time-stepping. Let us denote by

x n+1 p , ω n+1
p the locations and strengths of the particle at the end of the time step. Particles are then remeshed on the original grid points by interpolation, with the formula :

ω n+1 p = p ω n+1 p Λ x n p -x n+1 p h . ( 22
)
Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

In this study, the function Λ in the above formula is obtained by tensor product of a piecewise cubic third order interpolation kernel, the so-called M 4 kernel originally proposed in [START_REF] Monaghan | Extrapolating B splines for interpolation[END_REF] :

Λ(x) =              0 if |x| > 2 1 2 (2 -|x|) 2 (1 -|x|) if 1 ≤ |x| ≤ 2 1 -5x 2 2 + 3|x| 3 2 if |x| ≤ 1 (23)
The same kernel is used to interpolate grid velocity values onto particles in the time-stepping scheme to solve [START_REF] Nield | Convection in Porous Media[END_REF]. One of the main feature of particle methods is that they do not have a CFL restriction on the time step ∆t. It is indeed shown in [START_REF] Cottet | High order semi-lagrangian particles for transport equations : numerical analysis and implementation issues[END_REF] that a remeshed particle method for a transport equation with velocity u is stable under the condition

∆t ≤ C ∇u ∞ , ( 24
)
where C is a constant that only needs to be smaller than 1.

Force calculation

Assume a flow in R N around a bounded body S with velocity and vorticity vanishing at infinity and having a prescribed value u b on the body. The formula ( 15) is based on the velocity-pressure formulation of the penalization method. It assumes that the velocity is of the order of K inside the body. In a vorticity formulation however, there is no direct control of the velocity in the body. In practice one observes that the velocity in the body is small but larger than K. This can be seen in Figure 2 (left). In that case, we found that accurate computations of forces could not be achieved using the formula [START_REF] Rossinelli | Gpu accelerated simulations of bluff body flows using vortex particle methods[END_REF] and we instead used force evaluations in terms of vorticity distributions around the body.

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

The force F exerted by this flow on the body is given by the following formula :

F = - d dt S u dx, (25)
where S denotes the domain outside S. Integration by parts allows to express this force in terms of the vorticity, with N denoting the space dimension :

F = - 1 N -1 d dt S x × ω dx - 1 N -1 d dt ∂S x × (n × u)ds, (26)
where the first integral is the linear hydrodynamical impulse and the second integral is the bound vorticity (see [START_REF] Saffman | Vortex Dynamics[END_REF]). For an impulsively started obstacle at early times, the vorticity remains essentially attached to the body, and this formula can be used by taking a computational domain large enough to contain the vorticity support. When vorticity shedding occurs, this is no longer possible. Alternatively, one may use the Navier-Stokes equations and the expression [START_REF] Cottet | High order semi-lagrangian particles for transport equations : numerical analysis and implementation issues[END_REF] to obtain :

F = -ν ∂S ω × n ds + ∂S pn ds. (27)
In the case of a 2D circular cylinder of radius R, this yields the following formula [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] :

F = νR 2π 0 ω - → e θ dθ - 2π 0 ∂ω ∂n - → e θ dθ . (28)
This expression assumes that vorticity vanishes at infinity, which in practice is always the case, and allows to distinguish friction (in the first integral) and pressure drag (in the second integral).

For a detached bluff body flow this formula can be used as long as vorticity and its normal derivatives are accurately captured at the boundary. However, the penalization method just described is derived from a method which is only first order for the velocity. It therefore does not capture vorticity values and fluxes at the boundary with enough accuracy to give reliable force evaluation through [START_REF] Rasmussen | A multiresolution remeshed vortex-in-cell algorithm using patches[END_REF] with reasonable grid sizes. It is therefore necessary to use a method which does not involve integrals of the physical quantities on the volume of the body. The technique chosen here combines time derivatives, volumes integral of the vorticity momentum -like in formula (26) -and vorticity boundary values as well as fluxes on a surface away from the body, where these quantities are accurately computed. We follow the approach proposed in [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], and rewrite the force according

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls to the momentum equation as follows :

F = - d dt V u dx + Σ n • γ dx, ( 29
)
where V is a domain containing S, Σ denotes its boundary and

γ = ( 30
)
1 2 |u| 2 I - 1 N -1 u • (x × ω) + 1 N -1 ω • (x × u) - 1 N -1 x • ∂u ∂t I -x • ∂u ∂t + 1 N -1 [x • (∇ • T )I -x(∇ • T )] + T .
In the above equation T denotes the stress tensor :

T = µ(∇u + ∇u t ), (31)
where µ is the viscosity. For a sake of simplicity we develop this formula in the two-dimensional case where the computational box D and the volume control V are two nested rectangles. We use the notations given in Figure 1. We denote the velocity by u = (u, v) t . The streamwise and vertical Prepared using fldauth.cls components of the force can then be written computed as :

F x = - d dt V u(1 -χ S )
+ A1 A4 - 1 2 (v 2 -u 2 ) + uωy + y ∂v ∂t - 1 Re 2 ∂ 2 v ∂y 2 + ∂ 2 v ∂x 2 + ∂ 2 u ∂x∂y y - 2 Re ∂u ∂x dy,
and For a two dimensional viscous flow past a circular obstacle of diameter d with free stream velocity u ∞ , drag and lift coefficients are expressed by the following formulas :

F y = - d dt V v(1 -χ S ) dxdy + A2 A1 1 2 (v 2 -u 2 ) -vωx + x ∂u ∂t - 1 Re 2 ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 + ∂ 2 v ∂x∂y x - 2 Re ∂v ∂y dx + A3 A2 -uv + uωx + x ∂v ∂t - 1 Re 2 ∂ 2 v ∂y 2 + ∂ 2 v ∂x 2 + ∂ 2 u
C D = 2F x ρu 2 ∞ d , C L = 2F y ρu 2 ∞ d . ( 32
)

VALIDATION FOR FLOW PAST A TWO-DIMENSIONAL CIRCULAR CYLINDER

The vortex penalization has been validated in 2D and 3D fluid structure interaction problems at low Reynolds numbers [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF], in the context of multi-resolution [START_REF] Rasmussen | A multiresolution remeshed vortex-in-cell algorithm using patches[END_REF] and in GPU implementations for 2D flows [START_REF] Rossinelli | Gpu accelerated simulations of bluff body flows using vortex particle methods[END_REF]. To evaluate the accuracy of the method, in particular with respect to the value of the penalization parameter, we consider here validations on the classical benchmark of the flow around a two-dimensional circular cylinder for a wide range of Reynolds numbers.

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

There is an important literature about the classical benchmark of the 2D circular cylinder. This benchmark is covered by an important literature. A thorough survey of this subject can be found in [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF] [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF].

Vortex methods using vorticity flux boundary conditions have allowed to obtain reference results for a wide range of Reynolds numbers [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. Recently, the penalization-vortex method has been implemented on GPU processor. The accuracy and performance of this implementation against traditional CPU implementations has been tested on this benchmark in an explicit formulation of the penalization. We use here the same benchmark to illustrate the efficiency of the implicit formulation (equation ( 20)).

Numerical setup and flow conditions

The computational domain and the geometrical setup are shown in Figure 1, where D is a rectangle delimited by its boundaries Γ D . The non dimensional diameter d of the circular cylinder is equal to 1 and the body is centered at the origin. The whole computational domain is covered by a uniform Cartesian grid.

In this section, periodic boundary conditions on the box boundary Γ D were prescribed for the system [START_REF] Chorin | Vortex sheet approximation of boundary layers[END_REF]. A correction of the velocity field was performed at each time step in order to satisfy free stream velocity u ∞ = (u ∞ , 0) = (1, 0) imposed at the inlet. In all cases the size of the computational box was chosen large enough so that the effect of these boundary conditions was found to be negligible.

Convergence studies

First, we consider the case of an impulsively started cylinder flow with Re = 550 and study the convergence of the method with respect to the penalization parameter and the grid-size. The left picture of Figure 2 shows the mean velocity profiles inside the cylinder, averaged over the time

window T = [0, 3].
The results corresponding to the implicit penalization scheme for three values of λ are compared to the profile obtained with the explicit formulation for a mesh size h = 0.005.

The explicit formulation of ( 18) requires a maximum value of λ = 1/∆t to guarantee stability.

Copyright c 0000 John Wiley & Sons, Ltd. One can see from this experiment that, to obtain small values of the velocity inside the body, it is necessary to use larger value of λ, and therefore to implement the implicit scheme [START_REF] Bruneau | Passive control of the flow around a square cylinder using porous media[END_REF]. As already

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls mentioned, one can also notice that velocity values inside the body remain in all cases larger than 1/λ. Unlike for the velocity-pressure formulation, indeed, in the vorticity formulation there is no direct relationship between the velocity in the body and the parameter λ.

To have a more precise idea of the sensitivity of the results to the value of the penalization parameter, a λ-convergence study is carried out at Re = 550 with h = 0.005 and ∆t = 0.0025. The value of λ ranges from 1 to 10 6 , with λ = 10 8 as reference case (solid case). We evaluate the

L 1 , L 2
and L ∞ norms of the error e(y) in the mean velocity magnitude profile | u| (at x = 0) with respect to the reference case :

e(y) L 1 = ymax y min | u| λ=10 8 (y) -| u|(y) dy (33)
e(y)

L 2 = ymax y min | u| λ=10 8 (y) -| u|(y) 2 dy 1/2 (34) e(y) L ∞ = sup y∈[y min ,ymax] | u| λ=10 8 (y) -| u|(y) (35)
As Figure 3 (a) shows, the numerical error decreases as O(λ -1/2 ) when λ ≤ 10 2 and as O(λ -1 ) when λ ≥ 10 2 . The same study was performed at Re = 100 with a smaller grid size and time step, respectively equal to h = 0.0025 and ∆t = 0.001. The results are given by Figure 3 (b) and show that the error behaves as O(λ -1/2 ) when λ ≤ 10 4 and O(λ -1 ) when λ ≥ 10 4 . These convergence rates are in agreement with the numerical results found out in [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] where it is stated that to increase λ to high values (λ ≥ 10 6 ) a mesh refinement is necessary. Moreover, these numerical results are slightly better compared to the theoretical result given in [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], attesting that the upper bound on the global

L ∞ error of the penalization is u λ -u ≤ Cλ -1/2
. Furthermore, we can see on this figure that the circular cylinder reaches the solid state from λ = 10 4 , where e(y) ≤ 10 -4 . Thus, a permeable medium can be defined by setting the λ-values between 1 and 10 3 , with very low permeability for the latter one.

Finally, the right picture of Figure 2 shows a grid-refinement study performed at Re = 550 for a section located in the wake at x = 1.5 and a penalization parameter λ = 10 10 in the implicit formulation. It depicts the instantaneous velocity magnitude at time T = 3 for different grid sizes

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls and indicates that for this flow one can consider that convergence is reached for h ≤ 0.005.

Additional validations are given in the next section.

Flow analysis and validations

The first part of the study is related to the low Reynolds number analysis where the viscous effects are predominant. Following [START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodies of general geometry[END_REF] the adimensional time-step ∆t is determined by the condition ∆t/(h 2 Re) ∼ O(1). Various tests have been carried out by increasing the blockage ratio 1/L y , where L y denotes the heigth of the computational box, as shown on Figure 1. The size of the computational domain was chosen such that the effects of the boundaries on the shedding frequencies were negligible. Moreover, 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. This with 3250 × 1500 grid points. This corresponds to a mesh size h = 0.01. The penalization was introduced using the implicit formulation [START_REF] Bruneau | Passive control of the flow around a square cylinder using porous media[END_REF] and the penalization parameter was set to λ = 10 10 .

For Re < Re crit = 49, the flow regime is laminar and the solution is steady and stable (see [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF]).

For larger Reynolds numbers, the flow regime is still laminar but the steady solution becomes unstable (Re > Re crit ). On the left-hand side picture of the Figure 4 estimate of the mean relative error between the two curves yields

¯ = 1 N N |St ref -St| ≈ 0.6%,
where N denotes the number of points on each curve (i.e N = 13).

The evaluation of forces has been carried out using the momentum equation [START_REF] Russell | A cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow[END_REF]. The mean values and the amplitudes of the fluctuations of drag and lift coefficients for flows at Re = 100

and Re = 200 are computed and given in Table 1. These results are very close to experimental and numerical data collected by Russell & Jane Wang [START_REF] Russell | A cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow[END_REF] and Braza et al. [START_REF] Braza | Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder[END_REF].

To address a transitional case we focus now on an impulsively started flow at Re = 550 and compare our results to those of Koumoutsakos & Leonard [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] and Ploumhans & Winckelmans [START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodies of general geometry[END_REF]. The time evolution of the drag coefficient is reported in Figure 5 (a). In our computations the the grid convergence in a computational domain [-4, 8] × [-5, 5] is achieved with parameters h = 0.005, ∆t = 0.005 and λ = 10 8 . The various results are in excellent agreement. [START_REF] Apte | A numerical method for fully resolved simulation (frs) of rigid particleflow interactions in complex flows[END_REF][START_REF] Mittal | Flow-induced vibrations of a light circular cylinder at reynolds numbers 1000 to 10000[END_REF] for long time simulations. All these comparisons confirm an excellent agreement between our results and the reference results in the literature.

For Re = 3000 the flow is unstable and the shedding generates complex vortex pairings (see [START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF]). In this case, grid convergence was obtained in our method for λ = 10 8 , h = 0.0025 and ∆t = 0.0025. Figure 5 (c) shows again a very good agreement in the drag history between our results and those in [START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodies of general geometry[END_REF] and [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. In [START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodies of general geometry[END_REF] h is about 0.002. Vortex contours in Figure 6 (b) confirm this agreement.

The last simulation addresses the flow around the impulsively started cylinder at Re = 9500. Let us emphasize that a 2D validation study at such a regime makes sense since we only consider here the impulsive start of the flow, before the onset of 3D instabilities. In order to capture small spatial structures and the thin boundary layer, the grid size was reduced to the value h = 0.00125 in the same computational domain as the one used in the previous cases. The time step was set to 0.002 and λ = 10 8 . The comparison of our results to those obtained by Rasmussen et al. [START_REF] Rasmussen | A multiresolution remeshed vortex-in-cell algorithm using patches[END_REF] as well as Koumoutsakos & Leonard [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] is reported in Figure 5 (d). In [START_REF] Rasmussen | A multiresolution remeshed vortex-in-cell algorithm using patches[END_REF] a multilevel vortex method was used with a minimal grid resolution corresponding to h min = 1/1024. The discrepancy observed on Figure 5 (d) at early times of the drag history indicates that the penalization method has difficulties to establish the strong initial vortex sheet on the body of the cylinder. For later times our results are in excellent agreement with [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Rasmussen | A multiresolution remeshed vortex-in-cell algorithm using patches[END_REF].

To conclude this section, Table 3 is given to confirm that the choice of the grid resolution h at each flow regime is appropriate to capture the boundary layer thickness, expressed by δ ≈ D/ √ Re.

The validation study just performed for the 2D circular cylinder benchmark enables us to determine the appropriate grid resolutions and time steps needed to obtain accurate results in the more challenging cases that will be considered in the following.

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls The main problem encountered with the aerodynamics of a VAT architecture turbine is the complex unsteady flow phenomena occurring during its working cycle. For each revolution the blades undergo a highly unsteady relative motion with large variation of incidence and absolute value of the relative velocity. This complexity in the flow motion is mainly characterized by flow detachment, vortex shedding, dynamic stall and turbulence.

The study of the turbulent flow past vertical axis turbine has been addressed recently by some authors [START_REF] Shamsoddin | Large eddy simulation of vertical axis wind turbine wakes[END_REF], [START_REF] Tescione | Near wake flow analysis of a vertical axis wind turbine by stereoscopic particle image velocimetry[END_REF], [START_REF] Mcnaughton | Turbulence modelling of low reynolds number flow effects around a vertical axis turbine at a range of tip-speed ratios[END_REF] ; however the basic characteristics of the turbulent wake for such a kind of flow remains essentially an open issue. [START_REF] Cottet | A comparison of spectral and vortex methods in threedimensional incompressible flows[END_REF][START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high reynolds numbers[END_REF] provide validation studies of remeshed vortex methods against reference spectral methods in 3D laminar and turbulent flows. These studies give confidence that remeshed vortex methods, combined with penalization methods, can be efficient and flexible tools to handle VAT flows.

The VAT is rotating at a given velocity Ω and is immersed in a stream with a uniform far field velocity u ∞ . The blade radius and the airfoil chord are denoted by R and c, respectively (Figure 7 (a)).

The motion of a solid obstacle can be described with the vortex-penalization method as follows [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF]. S = S(t) is a moving multi-connected solid region. If one denotes by φ a level-set function which is negative inside S and positive outside, the characteristic function of S is defined by

χ = H(φ)
, where H(•) is a suitable step function. Moreover, in a rigid body motion, the level-set function has to satisfy

∂φ ∂t + (u s • ∇) φ = 0. ( 36
)
Copyright c 0000 John Wiley & Sons, Ltd. This bench study mainly aims to highlight the great capability of the present method to supply realistic results for moving bodies with unfitted boundaries and the grid resolution handled here is convenient to confirm the qualitative physical behaviour of computations.

Prepared using fldauth.cls c R Y X 0 Ω ∞ u c (a) (b)
For the Reynolds number Re = 1000, a wind turbine with a blade section chord c = 0.2 m, radius R = 0.3 m, rotating speed Ω = 0.185 rad/s and immersed in a stream with u ∞ = 0.025 m/s is considered (where ν air = 1.6 10 -5 m 2 /s at 20 Celsius degrees). For this Reynolds number, which corresponds to a low transitional flow with viscous effects, the time step is taken as ∆t = 0.001. A higher Reynolds number for a wind turbine operating at low-wind speed is studied by considering this VAT geometry with the working conditions Ω = 1.85 rad/s and u ∞ = 0.25 m/s.

The corresponding Reynolds number is Re = 10000 and the integration time step is chosen to be ∆t = 0.0002. On Figures 8 (a and Re = 10000. The isolines are displayed in the relative frame of reference, that is centered on the turbine axis and moves jointly with the blades. For the lower Reynolds number Re = 1000 the average is computed for a simulation time of ≈ 650s, which corresponds to ≈ 19 loops of the blade.

As the Figure 9 (a) shows, for this flow the recirculation zones (one for each blade) are concentrated between the rotation axis and the inner surface of the blade. For the higher Reynolds number Re = 10000, the mean flow is computed for a simulation time of ≈ 65s, which still corresponds to ≈ 19 cycles. In this case, for each blade three recirculating areas can be identified, as in Figure 9 Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls (b) : two zones are close to the airfoil (outer and inner surfaces) and the third region is the average in time of the wake generated from the trailing edge of each blade due to complex vortex interactions.

Finally, we study the time evolution of the enstrophy and the number of vortex particles in order to better analyze the properties of these two flow regimes. The enstrophy is a quantity that measures the dissipation of energy in terms of kinetic energy. In a 2D flow field the enstrophy is given by In the present cases the flow regime is characterized by unsteady separations generated by the periodical interaction of the wake with the blades. For Re = 1000, we notice that the principal vortex shedding frequency roughly corresponds to 2Ω, that is twice the angular velocity of the two-bladed VAT. For Re = 10000, different shedding frequencies can be distinguished from the rotational frequency of the blade. The investigation on the shedding frequencies is meaningful because it allows to relate the interaction of the vortex wake with the structure of the turbine. If some vortex shedding frequencies match the proper frequencies of the turbine, the structure can resonate with a self-sustained motion and the turbine is liable to aeroelastic phenomena. To evaluate the principal frequencies of vortex shedding we performed a discrete Fourier analysis of the enstrophy signal. In Figures 11 (a In Figure 11 (a) the shedding principal harmonic at 2 times the angular velocity of the turbine is confirmed : the wake interacts with a blade two times for each cycle of the turbine. Further couplings with the principal frequency can be identified at 4, 6, 8. Some high energy contributions are visible at frequencies lower than 2. A reason may be that the flow is not fully developed and did not achieve a statistical steady state yet.

Z = D |ω| 2 dσ
In Figure 11 (b) the spectrum exhibits shedding frequencies at 2 and 4 times the angular velocity.

The enstrophy signal is noisy and is characterized by random peaks at low frequency (Figure 10 (b)).

Following [START_REF] Pope | Turbulent Flows[END_REF], since the turbulent dissipation rate is equal to the turbulent kinetic energy production, the ratios of the Kolmogorov scale to largest scales are determined readily by means of the approximation ∼ U 3 L and the results are : The horizontal axis reports the discrete frequency ω i divided by the angular velocity Ω of the turbine.

η ∼ ν 3 L U 3 1 4 , τ η ∼ νL
Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

EXTENSION TO SOLID-POROUS-FLUID MEDIA AND FLOW CONTROL

This section is devoted to the study of flows in solid-porous-fluid media. It illustrates the efficiency of the vortex penalization method to handle problems where the physical properties or boundary conditions vary in the body.

Modeling and numerical setup

The penalization method has been used to model the porous-fluid-solid problems with primitive pressure-velocity formulations [START_REF] Bruneau | Numerical modelling and passive flow control using porous media[END_REF][START_REF] Bruneau | Passive control of the flow around a square cylinder using porous media[END_REF][START_REF] Bruneau | Control of vortex shedding around a pipe section using a porous sheat[END_REF]. This work is dedicated to the extension of this approach to the vorticity equations and to different configurations. We consider a solid body surrounded by a porous layer and immersed in an incompressible fluid. In such a system, on can identify five in Figure 12 (left)). The goal is thus to compute the flow simultaneously in the porous medium, the fluid region and at the interface between the two media. Several approaches have already been proposed to handle this problem. The first one, which is mainly used to model flows past a solid with permeable walls [START_REF] Hahn | Direct numerical simulation of turbulent channel flow with permeable walls[END_REF][START_REF] Jimenez | Turbulent shear flow over active and passive porous surfaces[END_REF], avoids to solve the porous flow by enforcing appropriate porous-fluid boundary conditions [START_REF] Beavers | Boundary conditons at a naturally permeable wall[END_REF]. This method neglects the physics of the porous medium and does not allow to have an overall view of the problem. A second approach considers that the porous flow has to be accurately computed, and solves the governing equations in each region by coupling the Darcy equations and the Navier-Stokes equations through an appropriate condition at the interface [START_REF] Breugem | The laminar boundary layer over a permeable wall[END_REF][START_REF] Hanspal | Numerical analysis of coupled stokes/darcy flows in industrial filtrations[END_REF]. This approach is very delicate to handle especially because of the requirement

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls of this appropriate boundary condition. The method presented here, based on the vortex-penalization technique, considerably simplifies the task since a unique equation is used for the whole domain and enables to accurately model each of the different regions by varying the value of the dimensionless penalization factor λ. In 2 dimensions, the penalized vorticity equation reads :

∂ω ∂t + (u.∇)ω = 1 Re ∆ω + ∇ × [λ χ S (u s -u)]. (37)
We recall that λ essentially depends, in the inverse proportion, on the physical intrinsic permeability of the medium k which directly derives from the Darcy law (equation ( 7)). Varying the value of λ thus directly characterizes the different media.

Modeling the physics of three different regions enables one to deal with engineering problems involving porous media. In the following, this approach is validated for a simple but informative passive flow control problem : the solid-porous-fluid configuration is applied to cover a semi-circular cylinder geometry with a porous coating. The porous coating is applied on the obstacle external surface in order to modify the vorticity generation of the boundary layer and the vortex shedding.

The presence of a porous medium at the solid-fluid interface indeed imposes a kind of mixed boundary condition intermediate between the no-slip and the slip one on the solid boundary [START_REF] Carbou | Brinkmann model and double penalization method for the flow around a porous thin layer[END_REF]. As a result, the shear forces are reduced and the flow dynamics is smoothed [START_REF] Bruneau | Passive control of the flow around a square cylinder using porous media[END_REF][START_REF] Bruneau | Control of vortex shedding around a pipe section using a porous sheat[END_REF][START_REF] Bruneau | Numerical modelling and passive flow control using porous media[END_REF][START_REF] Bruneau | Passive control around a two-dimensional square back ahmed body using porous media[END_REF]. This control technique, which allows to keep the obstacle geometry unchanged, is expected to reduce drag forces and vortex induced vibrations, thereby improving the aerodynamic properties of the obstacle. The semi-circular cylinder can be considered as a simplified section of an outside rear-view mirror of a vehicle. The mirrors, due to their spanwise location, indeed generate a non-negligible wake which interferes with the vehicle sides. This provides a good motivation to perform flow control past such obstacles. Moreover, as it was shown in [START_REF] Rouméas | Analysis and control of the near-wake flow over a square-back geometry[END_REF][START_REF] Bruneau | Coupling active and passive techniques to control the flow past the square back ahmed body[END_REF], a flow past a square back obstacle is not dominated by longitudinal vortical structures, therefore a preliminary two-dimensional study can be useful to supply information on general trends for a control.

In the following, the semi-circular cylinder has a total dimensionless diameter of Figure 12 (right)). In the sequel we will consider flows both at transitional (Re = 550) and highly transitional (Re = 3000) regimes.

Validation and convergence studies for flow past a solid semi-circular cylinder

Compared to the circular cylinder considered in the previous section, the semi-circular cylinder presents singular points which make its simulation more challenging. In order to verify the accuracy of the method in the case of flow past a solid semi-circular cylinder, we proceed like in section 3.

We first compare the results obtained using the present method with the ones of the literature, for different Reynolds numbers. Secondly we present the convergence of the method with respect to the penalization parameter and the grid-size.

Validation

The benchmark of incompressible flow past a semi-circular cylinder is not a classical one and only few references already exist on this subject. Nevertheless, we aim at least to compare our results to numerical and experimental ones. To do so, we respectively consider the works of Farhadi et al. [START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF] and Boisaubert and Texier [START_REF] Boisaubert | Effect of a Splitter Plate on the Near-Wake Development of a Semi-Circular Cylinder[END_REF]. The flow analysis presented hereafter are based on the near-wake

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls (green triangles) and experimental data [START_REF] Boisaubert | Effect of a Splitter Plate on the Near-Wake Development of a Semi-Circular Cylinder[END_REF] (black squares).

length measurement, the Strouhal number and the drag and lift coefficients. The Reynolds number ranges between 100 and 300 in this study.

As can be seen in Table 4, the results obtained at Re = 100 and Re = 300 show a good agreement with the numerical values given by [START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF]. A more detailed comparison is performed at Re = 200

and the related results are reported in Figure 13. As regards drag and lift coefficients, one can say according to Figure 13 (a) that our results coincide very well with those of [START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF] since we achieve cD = 1.85 ± 0.15, cL = 0 ± 0.53 and S t = 0.230 against cD = 1.85 ± 0.12, cL = 0 ± 0.50 and S t = 0.215 in [START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF]. Figure 13 (b) gives a comparison of the normalized near-wake length L with the numerical results of [START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF] and the experimental data of [START_REF] Boisaubert | Effect of a Splitter Plate on the Near-Wake Development of a Semi-Circular Cylinder[END_REF]. The near-wake length is defined as the x-coordinate of the point located on the centreline in the x-direction where the streamwise component of the velocity vanishes. As can be seen on the figure, the values obtained with the present method greatly coincide with the numerical and experimental reference results.

Convergence studies

In this section we present convergence studies of the vortex-penalization method for the flow past

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls a solid semi-circular cylinder at Re = 550 and Re = 3000. We first focus on the convergence of the present method with respect to the penalization parameter λ. The way this convergence study is carried out is exactly the same as presented in section 3.2. As Figure 14 shows, the semi-circular cylinder reaches the solid state from λ = 10 5 , where e(y) ≤ 10 -4 , for both regime. A permeable area can thus be defined with λ-values lower than 10 4 . In the oncoming control study, we will vary the penalization parameter until 10 3 inside the added porous layer in order to model a permeable medium.

We now turn to a grid convergence study. It is performed on three grid levels. For Re = 550

we consider three mesh sizes are : h = 0.01, h = 0.005 and h = 0.0025 in the computational [START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF][START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF]. For Re = 3000 as the boundary layer is thiner, we perform the convergence study with three finer mesh sizes, respectively h = 0.005, h = 0.0025 and h = 0.00125.

domain D = [-4, 8] × [-
The simulations are carried out in a smaller domain D = [-2.5, 5] × [-3, 3] in order to reduce the computational cost induced by the finest grid. The mean drag and enstrophy values are reported in Table 5. At Re = 550, on the basis of these results, one can consider that the grid convergence is achieved for h = 0.005 with ∆t = 0.0025. At Re = 3000 grid convergence is obtained with h = 0.0025 and h = 0.00125, which leads us to choose h = 0.0025 and ∆t = 0.002 for the oncoming simulations. Let us emphasize the interest of the vortex methods which enables here to handle a highly transitional regime using a relatively large time step compared to classical Eulerian methods.

Parametric study with respect to the penalization parameter λ

In this section we study the influence of the porous layer permeability on the flow behavior and the efficiency of such a passive control. This parametric study is performed at Re = 550 and Re = 3000

with different values of the porous permeability inside the layer, ranging from λ = 0.125 (very high permeability) to 10 3 (very low permeability). The results are compared to the one corresponding to the solid (λ = 10 8 ) and fluid (λ = 0) cases. In a first place, the thickness of the coating is set to ten percent of the diameter : τ = 0.1 (see Figure 15) for all the cases. According to the grid convergence performed in the previous section (Table 5), the numerical simulations are performed

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls (see Figure 12 (right)).

In order to analyze the effects of our control approach we compare global flow quantities like the drag force (F x ), computed according to the momentum equation [START_REF] Russell | A cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow[END_REF], the enstrophy (Z) allowing to measure the dissipation effects in the flow as well as the delay of transition to turbulence, and the Strouhal number (S t ). Note that in this study we consider drag forces instead of drag coefficients since the C D formula (32) involves the diameter d of the obstacle which is not clearly defined because of the porous coating.

The results of the parametric study performed at Re = 550 are given in Table 6. Concerning the Strouhal number S t , one can notice that it continuously increases as the λ value decreases. The presence of a permeable layer around the obstacle indeed enables to raise the shedding frequency, and in the limit λ → 0 one tends to the flow dynamics of the fluid case which is characterized by smaller vortices swirling at higher frequency (Figure 17 (a)). Regarding Fx , Table 6 shows that, as expected, the drag reduces as the controlled case tends to the fluid case. Nevertheless, one can In fact, this range of permeability values is responsible for a high resistance of the flow inside the layer inducing higher friction forces at the interface, and thus a higher total drag. The quantitative results reported in Table 6 are also highlighted in Figures 18 and 20 respectively depicting the mean velocity and mean vorticity fields. As already emphasized in the literature [START_REF] Bruneau | Coupling active and passive techniques to control the flow past the square back ahmed body[END_REF][START_REF] Bruneau | Passive control around a two-dimensional square back ahmed body using porous media[END_REF], the main source of drag forces for bluff bodies with square back is the low pressure recirculation zone in the near wake of the vertical wall. With λ = 1, the near wake structures are smaller and the back recirculation zone is drastically reduced (Figure 20). Moreover, the mean velocity of the wake clearly diminishes (Figure 18) and the structure of the global and near wake resembles the one of the fluid case. This feature explains the important drag reduction observed in Table 6 for λ = 1 (-27%) and in general for small values of λ. On the contrary, the flow behavior obtained with λ = 10 looks very similar to the one of the uncontrolled case with an important transversal extent of the wake (Figure 20)

and high velocities in the flow compared to the uncontrolled case (Figure 18). This confirms the results given in Table 6 and allows us to conclude that this value of λ is not suitable at this regime to achieve drag reduction. However, with λ = 10 we obtain an important reduction of the enstrophy which implies an efficient wake regularization and thus a delay in flow transition. This observation leads us to the analysis of control efficiency in terms of enstrophy, for which we observe a more complex behavior. Indeed, as can be read in Table 6, the enstrophy decreases with λ to reach an optimal value, for λ = 2.5 (-41% compared to uncontrolled case). Below this λ value, the enstrophy

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls increases again to finally reach the value corresponding the reduced obstacle (fluid case). Therefore, when λ ≤ 10, then Zλ < Zfluid , which demonstrates the non-monotonic behavior of the enstrophy, in the range of permeabilities taken into account in this parametric study (for all values of λ in {0.125, 1000}, Zλ < Zsolid ). Figure 16 (a) summarizes these observations made from Table 6 and shows the evolution of the mean drag force and the mean enstrophy with respect to the λ parameter.

At Re = 3000, the evolution of the Strouhal number as a function of λ is similar to the one observed at Re = 550 : the more permeable the added layer is, the higher the shedding frequency is. This behavior is represented by the vorticity fields of Figure 17 (b), where one can clearly see the benefits of the porous coating which enables to go from a chaotic wake with large eddies of strong intensity (uncontrolled case) to an axisymmetric Von Kármán street composed of numerous small vortices, similar to the one of the fluid case (controlled case with λ = 1). On the other hand, the results obtained in terms of mean drag force and mean enstrophy are also comparable to those obtained at Re = 550, with a more pronounced behavior. Concerning the drag force, we indeed notice on Table 7 and Figure 16 study, we identify an optimum λ value, still equal to 2.5, for which the enstrophy reaches a global minimum value leading to a reduction of 44% (Table 7 and Figure 16 (b)). In a general way, this nonmonotonic behavior of the enstrophy proves that the regularization effects reported here are not due to a change of the Reynolds number, since we have for instance Zλ=2.5 < Zλ=1 < Zfluid , whereas Re λ=2.5 > Re λ=1 > Re fluid . In fact, the introduction of a porous layer generates Kelvin-Helmoltz instabilities at the solid-porous interface. These instabilities, depicted on Figure 22 at Re = 3000

for the λ = 1 test case for example, change the shear forces at the interface and therefore modify the vortex shedding mechanism and the global flow dynamics. The presence of such instabilities explains in particular the qualitative and quantitative differences between the fluid case and the very

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls permeable cases.

Finally, the most important observation which emerges from this parametric study is the difference of behavior between the enstrophy and the drag evolution. We indeed highlighted the systematic enstrophy reduction in the presence of a porous coating, regardless of its permeability, whereas some λ values were found to be responsible for a non negligible drag increase, especially at Re = 3000. This is directly due to the way these two quantities are evaluated. The enstrophy quantifies the energy dissipation of the whole system, it is therefore evaluated in the entire domain.

Thus, the laminar velocities caused by the presence of the permeable layer and appearing in the vicinity of the detachment points, lead to a delay in transition to turbulence and directly impact the transport of the vorticity which is regularized. On the contrary, the drag force is evaluated on the walls of the obstacle or in the vicinity of the body (see section 2.3) and is directly linked to the vorticity generation. Therefore, the high friction phenomenon caused by some permeabilities (namely intermediate permeabilities) and implying a resistance of the flow in the front part of the obstacle, is directly reported in the drag diagnostic and in the very near-wake analysis. The λ = 10 test case at Re = 3000 clearly illustrates these statements : one can indeed notice the beneficial effects of the permeable coating on the wake, which becomes axisymmetric (Figure 17 (b) and does not show high velocity gradients anymore (Figure 19). As a consequence, we achieve a reduction of enstrophy (Table 7). On the opposite, the effects on the very near-wake and on the drag in particular are clearly not beneficial (Table 7) as can be seen on Figure 21. The case λ = 10 shows important recirculation zones in the rear back of the solid, responsible for a pressure decrease and thus for a drag growth. Therefore, the drag forces are more sensitive to the permeability values.

To conclude this λ-parametric study, we can assert that the control should be implemented with high permeability coefficients (λ ≤ 2.5) in order to combine both a significant drag reduction and a regularization of the wake.

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

Influence of the layer thickness τ on the flow control efficiency

In the previous section, the effects of the control were studied for different permeabilities using a constant layer thickness τ = 0.1. In this section, the effect of the thickness are investigated, and a special focus is made on thin layers in order to reduce the use of porous materials and to find out a compromise between practical engineering constraints and control efficiency. The same cases (Fig. 15) are now considered with the thickness of the layer set to τ = 5%d = 0.05 and to τ = 2.5%d = 0.025. Table 8 purpose of a companion paper [START_REF] Mimeau | Passive Flow Control Around a Semi-Circular Cylinder Using Porous Coatings[END_REF], where it is in particular shown that an efficient solution consists in locating the permeable zones only on top and bottom of the body, in the vicinity of the flow detachment points.

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

CONCLUSION

In this paper a method combining penalization, vorticity formulation and particle discretization was described and applied to complex fluid-solid configurations. The penalization approach greatly simplifies the implementation of vortex particle methods to handle the no-slip boundary conditions.

FIGURE 1 .

 1 FIGURE 1. Computational domain D and control volume V .

FIGURE 2 .

 2 FIGURE 2. Left picture : mean velocity magnitude profiles in a section across the center of the cylinder, at x = 0, obtained with implicit penalization (λ = 10 3 , 10 6 , 10 9 ) and explicit penalization (λ = 1/∆t) at Re = 550. Right picture : grid refinement study for the axial velocity at T = 3 on a cross-section at x = 1.5 behind the cylinder.

FIGURE 3 .

 3 FIGURE 3. λ-convergence study for flow past a circular cylinder at Re = 550 with h = 0.005 and ∆t = 0.0025 (a), and at Re = 100 with h = 0.0025 and ∆t = 0.001 (b). The errors e L ∞ , e L 1 and e L 2 are plotted against λ values.

  study shows a first-order convergence for the lift coefficient and Strouhal number and thirdorder convergence for the drag coefficient. It demonstrates that a domain extending to about 6 cylinder diameters downstream is sufficient to provide reliable flow diagnostics. The subsequent simulations have been performed in a computational box with dimensions [-7.5, 25] ×[-7.5, 7.5]

  , the amplitude A ω of the fluctuation of the vorticity on a monitoring point P = (2.5, 0.5) is plotted for various Reynolds numbers. For a flow regime close to the bifurcation point Re crit (Re < 60) the wake instabilities grow very slowly, so the oscillation study is started at Re = 60. Nevertheless, the exact critical Reynolds number is obtained with a linear extrapolation of the amplification factor curve. The righthand side picture of the Figure4shows the Strouhal number (given by the shedding frequency) as a function of the Reynolds number, obtained in our calculations and in the exprimental work[START_REF] Williamson | Vortex dynamics in the cylinder wake[END_REF]. AnCopyright c 0000 John Wiley & Sons, Ltd.Prepared using fldauth.cls

FIGURE 4 .

 4 FIGURE 4. Amplification factor Aω curve (left). Strouhal-Reynolds curve (right).

  Figure 6 (a) depicts the vorticity isocontours which are very similar to the vorticity field presented by Ploumhans & Winckelmans[START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodies of general geometry[END_REF] where h = 6.03 × 10 -3 .The next simulation was performed at Re = 1000. At this regime, advection effects become predominant in the flow field. The boundary layer thickness decreases, requiring a finer grid. The mesh size is consequently set to h = 0.0025 in a computational domain[-4, 8] × [-5, 5] and the time step is set to ∆t = 0.003. The non-dimensional time evolution of the drag coefficient is given in Figure5(b) and compared to the result in[START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF].

FIGURE 5 .

 5 FIGURE 5. Drag evolution for Re = 550 (a), Re = 1000 (b), Re = 3000 (c) and Re = 9500 (d). Comparisons of present method (solid red line) with Koumoutsakos & Leonard [2] (dashed line). Solid black line are results from Ploumhans & Winckelmans [8] in (a), (b), (c) and from Rasmussen et al. [28] in (d).

FIGURE 6 .

 6 FIGURE 6. Comparison of vorticity isocontours for the impulsively started cylinder between present method (left columns) and Ploumhans & Winckelmans [8] (right columns) at (a) T = 1, 3, 5 at Re = 550 and (b) T = 1, 2, 3 at Re = 3000.

FIGURE 7 .

 7 FIGURE 7. (a) Geometrical setup of VAT. (b) Distance function φ iso-contours ; the dashed rectangle Σ φ is the support of φ.

  ) and 8 (b) instantaneous vorticity isolines for different rotating times of the blades at Re = 1000 and Re = 10000 are respectively plotted. The varying incidence of the velocity occurring at each blade generates separation and vortex shedding ; as the figures show, the flow is highly unsteady and the blade strongly interacts with the wake. The flow for Re = 1000 generates an almost regular Von-Kármán-like wake with large structures. By contrast, the vortex shedding for the Re = 10000 becomes very quickly unsteady where small high frequency vortical structures interact with each other before being transported downstream.A complementary description of the flow around the turbine is given by the averaged-in-time streamlines for the absolute flow field shown on Figures9 (a) and 9 (b), for a wake at Re = 1000

  , where D is the computational domain. In figures 10 (a) and10 (b) the comparison between the time history of the enstrophy and the number of particles is shown for Re = 1000 and Re = 10000 respectively. As the Figure10 (a)shows, the amplitude of the oscillations remains quite small compared to the Re = 10000 regime (Figure10 (b)), where the frequency of the shedding appears through large oscillations. The figures show the high correlation between the time evolution of the enstrophy and the number of particles. In addition one can observe a characteristic frequency, roughly corresponding to 10 rotations of the blades, for which large vortices are shed outside the computational domain. Further computations would be needed to confirm and understand this phenomenon.

  ) and 11 (b) the spectra of the enstrophy versus the frequency of the signal are presented, where A enstrophy is the amplitude of the enstrophy and ω i /Ω is the nondimensional discrete frequency of the signal. The spectra does not report the mean value corresponding to the first Copyright c 0000 John Wiley & Sons, Ltd.Prepared using fldauth.cls coefficient i = 0 of the Fourier series. To avoid the effect of the initial transient stage, the Fourier analysis was applied on the signals after the first three cycles of the turbine.

U 3 1 2FIGURE 8 .FIGURE 9 .

 189 FIGURE 8. Iso-vorticity lines of the wake past the VAT at Re = 1000 (a) and at Re = 10000 (b). Snapshots respectively at times T = 3.8s, 11.4s, 18.9s, 26.5s and T = 0.38s, 1.14s, 1.89s, 2.65s.

FIGURE 10 .

 10 FIGURE 10. Time history of the dimensionless enstrophy (blue) and the number of vortex particles (red) for Re = 1000 (a) and Re = 10000 (b). The axis of the abscissas represents the number of the cycles performed by the turbine.

FIGURE 11 .

 11 FIGURE 11. Frequency spectrum of the enstrophy signal computed for Re = 1000 (a) and Re = 10000 (b).

  different flow regimes corresponding to different regions in the fluid. A first region is the boundary layer inside the porous medium, close to the solid wall. This region is very thin compared to the second region which is characterized by the homogeneous porous flow with Darcy velocity (regions numbered 1 and 2 in Figure 12 (left)). In the vicinity of the porous-fluid interface, two transient layers can be recognized (numbered 3 and 4 in Figure 12 (left)). The first one corresponds to the porous layer where the velocity increases to reach a value u i at the interface, and the second one to the fluid boundary layer leading to the free external flow. The size of the fluid boundary layer is determined by u 0 -u i where u 0 denotes the velocity of the main fluid flow (numbered 5

FIGURE 12 .

 12 FIGURE 12. Velocity profile in the vicinity of a porous medium (left). Computational domain for flow control past semi-circular cylinder (right)

FIGURE 13 .

 13 FIGURE 13. (a) Drag and lift evolution for flow past a semi-circular cylinder at Re = 200. (b) Normalized near-wake length L/d of a semi-circular cylinder at Re = 200 : comparison with numerical results[START_REF] Farhadi | Effect of a Splitter Plate on Flow Over a Semi-Circular Cylinder[END_REF]

FIGURE 14 .

 14 FIGURE 14. λ-convergence study for flow past a semi-circular cylinder : the errors e L ∞ , e L 1 and e L 2 are plotted against λ values at Re = 550 (a) and Re = 3000 (b).

FIGURE 15 .

 15 FIGURE 15. (left) Solid/uncontrolled case, (center) Porous case, (right) Fluid case.

  (b) that intermediate permeabilities (λ ≈ 10) completely deteriorate the performances. On the contrary, below λ = 5, the drag reduction is enhanced. This assessment is confirmed by Figures19 and 21showing in a case of a very permeable layer (λ = 1) a symmetric, thin and regular wake as well as low velocity gradients and low vorticity values in the vicinity of the detachment points and the recirculation zone. Regarding enstrophy, similarly as in the Re = 550

7 FIGURE 16 .

 716 FIGURE 16. Effects of various layer permeabilities on the mean drag force Fx (red curves) and the mean enstrophy Z (blue curves) for flow past a semi-circular cylinder at Re = 550 (a) and Re = 3000 (b) with τ = 0.1.

  reports the control performances in terms of global flow quantities, obtained for two controlled cases, λ = 10 and λ = 1, with τ = 0.05 and τ = 0.025. Generally, as expected, for both regimes (Re = 550 and Re = 3000) a slight decrease of control efficiency is observed when τ diminishes. Nevertheless, at Re = 3000, one can observe a change of trend. We indeed achieve a drag increase for both controlled cases (λ = 1, λ = 10), even if the enstrophy reduction is still maintained. Thus, to achieve a significant drag reduction it becomes necessary to think about the position of the permeable areas around the obstacle in order to reduce the high resistance currently experienced by the flow only in the front part of the layer. This question is the Copyright c 0000 John Wiley & Sons, Ltd.

FIGURE 17 .

 17 FIGURE 17. Instantaneous vorticity fields for flow past a semi-circular cylinder at Re = 550 (a) and Re = 3000 (b). The vortex shedding of the two controlled cases (λ = 10 and λ = 1, with τ = 0.1) are compared to the one of the uncontrolled case and the fluid case (reduced obstacle).

5 FIGURE 20 .FIGURE 22 .

 52022 FIGURE 20. Zoom of the mean vorticity fields and isolines for the flow past a semi-circular cylinder at Re = 550.

  It enables to consider the governing equations in the whole computational domain and enables the use of simple Cartesian grids and fast Poisson solvers to compute velocity fields. After a detailed validation on the circular cylinder benchmark, a the first application to moderate and high Reynolds flows around a moving vertical axis turbine was presented. The results showed the efficiency of this method to realistically model flows past moving obstacles. The second application illustrated the capability of the method to handle solid-fluid-porous media. The distinction between the three different media is indeed performed by varying the value of the penalization parameter without prescribing a boundary condition at the solid boundary or at the porous-fluid interface. An application to passive flow control past a semi-circular cylinder was performed at transitional and highly transitional regimes, consisting in adding a porous sheath on the obstacle surface in order to smooth the flow dynamics. The presence of a porous layer at the solid-fluid interface is responsible for a decrease of the shear forces and the vorticity generation of the boundary layer, leading to significant wake stabilization and drag reduction.

Table 2

 2

shows the comparison of mean drag Copyright c 0000 John Wiley & Sons, Ltd. Prepared using fldauth.cls coefficient, amplitude fluctuation of lift coefficient and Strouhal number with other studies

TABLE I .

 I Comparison of drag coefficient, lift coefficient, and Strouhal number for flow past a circular cylinder at Re = 100 and Re = 200.

			Re=100			Re=200
	Authors	cD	cL	St	cD	cL	St
	Braza et al. [30]	1.36 ± 0.015 0 ± 0.25	0.160	1.40 ± 0.050 0 ± 0.75 0.190
	Russel and Wang [29] 1.43 ± 0.009 0 ± 0.322 0.172	1.45 ± 0.036 0 ± 0.63 0.201
	Present method	1.40 ± 0.010 0 ± 0.32	0.165	1.44 ± 0.05	0 ± 0.75 0.200

TABLE II .

 II Comparison of mean drag coefficient, mean lift coefficient, and Strouhal number for flow past a circular cylinder at Re = 1000.

			Re=1000
	Authors	cD	cL	St
	Apte et al. [31]	1.50	0 ± 1.36 0.238
	Mittal et al. [32]	1.48 ± 0.21 0 ± 1.65 0.250
	Present method	1.51 ± 0.23 0 ± 1.54 0.245

TABLE III .

 III Justification of the choice of mesh size h based on the boundary layer thickness δ, for the circular cylinder validation study.

	Re	Boundary layer thickness (δ) Chosen mesh size (h) Number of mesh cells in boundary layer
	550	∼ 0.043	0.005	8.6
	1000	∼ 0.032	0.0025	12.8
	3000	∼ 0.018	0.0025	7.2
	9500	∼ 0.010	0.00125	8

TABLE IV .

 IV Comparison of mean drag coefficient and Strouhal number for flow past a semi-circular cylinder at Re = 100 and Re = 300.

			Re = 100	Re = 300
	Authors	cD	St	cD	St
	Farhadi et al. [50]	1.73 0.180	1.955 0.22
	Present method	1.70 0.192	1.90	0.24

Copyright c 0000 John Wiley & Sons, Ltd.

Prepared using fldauth.cls

C. MIMEAU, F. GALLIZIO, G.-H. COTTET, I. MORTAZAVI