Chlo É Mimeau
email: chloe.mimeau@imag.fr

Iraj Mortazavi
email: iraj.mortazavi@math.u-bordeaux1.fr

Georges-Henri Cottet
email: georges-henri.cottet@imag.fr

Passive Flow Control Around a Semi-Circular Cylinder Using Porous Coatings

In this paper, the passive control of flow past a semi-circular cylinder is investigated. This control is achieved by adding a porous medium between the solid obstacle and the incompressible fluid in order to reduce drag forces and regularize the flow. A vortex-penalization method is chosen to easily model the flow in the different media. Several configurations of the porous layer are investigated and parametric studies are performed in order to determine the most efficient passive flow control devices. This control study can be considered as a first step to propose efficient strategies to regularize the flow around a side-view mirror.

Nomenclature

C D drag coefficient D computational domain F, S fluid domain and solid domain F D drag force Re Reynolds number Z enstrophy d non-dimensional diameter h reference mesh size k intrinsic permeability l re f height of the obstacle u re f reference velocity u = (u, v) velocity field u s body rigid motion ū mean velocity magnitude u ∞ free stream velocity Γ D computational domain boundaries ∆t time step Φ porosity λ penalization parameter µ dynamic viscosity ν kinematic viscosity ρ density of the fluid τ porous layer thickness χ S characteristic function ω vorticity field

Introduction

On a ground vehicle, the outside mirrors, due to their spanwise position, indeed generate a non-negligible wake which interferes with the flow past car sides. They are responsible of up to 10% of the total vehicle drag but they only represent 0.5% of the total projected surface, which accounts for a good motivation to perform flow control past these obstacles. This work is devoted to the control of flow past a two-dimensional semi-circular cylinder which can be considered as a simplified section of a side-view mirror. 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 three-dimensional vortical structures, therefore a preliminary two-dimensional study can be useful to supply information and general trends for a further control study in three dimensions around a hemisphere. The aim is to use a control device easy to set up, low cost and allowing to keep the geometry unchanged. As active control devices can be hardly implemented in such a case, an efficient passive strategy seems to be a good alternative. A suitable solution has already been proposed by Bruneau and Mortazavi in [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 | Flow Manipulation Around the Ahmed Body With a Rear Window Using Passive Strategies[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]. It consists in adding a porous sheath on the obstacle surface in order to reduce the vorticity generation of the boundary layer. 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 decreased and the flow dynamics is smoothed. Consequently, the problem we have to solve involves three different media, namely the solid obstacle, the porous layer and the fluid. An easy way to tackle it is to use the penalization method [START_REF] Angot | A Penalization Method to Take Into Account Obstacles in Incompressible Viscous Flows[END_REF]. This method is based on a unique model, the Brinkman-Navier-Stokes equations, which are obtained by adding in the Navier-Stokes equations a penalization term, depending on the intrinsic permeability. Three values of this coefficient represent the three different media. This method can be easily implemented since it enables to consider the governing equations in the whole computational domain. Moreover it does not require to prescribe a boundary condition at the solid boundary or a condition at the porous-fluid interface. The penalization method has been recently implemented in vortex methods to deal with fluid-structure interaction problems [START_REF] Coquerelle | A Vortex Level Set Method for the Two-Way Coupling of an Incompressible Fluid with Collinding Rigid Bodies[END_REF][START_REF] Rossinelli | Gpu accelerated simulations of bluff body flows using vortex particle methods[END_REF][START_REF] Cottet | A Vortex Penalization Method for Flows with Moving Immersed Obstacles[END_REF]. It was extended to deforming bodies and applied to find optimal swimming modes of fishes [START_REF] Gazzola | C-start: optimal start of larval fish[END_REF][START_REF] Van Rees | Optimal shapes for anguilliform swimmers at intermediate reynolds numbers[END_REF]. In the present work, we consider a two-dimensional viscous and incompressible flow past a semi-circular cylinder. We also introduce some porous coatings on the obstacle in order to manipulate the flow. The Brinkman-Vorticity Transport Equations are solved using remeshed vortex methods [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF][START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Cottet | A Level-Set Formulation of Immersed Boundary Methods for Fluid-Structure Interaction Problems[END_REF] which are characterized by their efficiency and robustness at high Reynolds numbers since they allow to focus the computational task on the rotational zones and to solve the flow equations in a fast Lagrangian way.

In the following we first describe the methodology and show some validations and convergence studies. Then we present the results of two-dimensional passive flow control for several devices differing from each other in terms of the permeability, the thickness and the geometrical configuration of the added porous layer. A careful comparison is carried out between controlled and uncontrolled flows at transitional (Re = 550) and highly transitional regime (Re = 3000).

Vorticity Formulation And Vortex Methods

The dynamics of an incompressible flow 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, u the velocity and Re the Reynolds number. In this work, flow simulations are based on particle methods. The fluid particles which are displaced by convection and diffusion are characterized by their position and their vorticity. The vorticity transport is expressed by the Helmholtz equation (or Vorticity Transport Equation), obtained taking the curl of the incompressible Navier-Stokes equations (Eqn. 2) and given in 2D by :

∂ω ∂t + u.∇ω = 1 Re ∆ω in D, (3)
where ω denotes the vorticity. The Poisson equation

∇ 2 u = -∇ × ω, (4)
obtained from continuity equation (Eqn. 1), enables to recover velocity field once the vorticity field is known. The previous equations are approximated using a Vortex method [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF][START_REF] Chorin | Numerical Study of Slighly Viscous Flow[END_REF]. These methods are very robust and low-cost to simulate high Reynolds number recirculating flows since the computationel task is mainly focused on vortical structures (see for example [START_REF] Ghoniem | Vortex Simulation of Laminar Recirculating Flows[END_REF][START_REF] Sethian | Validation Study of Vortex Methods[END_REF]). In this kind of approach, the VTE equation (Eqn. 3) is solved using a two-fractional step (or viscous splitting) method. It relies on approximating separately the diffusion and convection terms at each time step. The convective part is solved using a "Vortex-In-Cell (VIC)" method (see e.g. [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF] and [START_REF] Creusé | Vortex Simulation of Active Control Strategies for Transitional Backward-Facing Step Flows[END_REF] for its application to active flow control) with a semi-Lagrangian resolution. In this fractional step a convective velocity is associated to each finite vortex element through a high order interpolation procedure, and the displacement is achieved using a Runge-Kutta method. In order to avoid Lagrangian distortion, particles are then remeshed on the original grid using the same interpolation kernel as the one used previously to interpolate grid velocity values onto the particles. Finally, for computational efficiency and accuracy, diffusion and Poisson equation are solved on the grid using Fast Fourier Transforms (FFT).

3 Vortex Penalization Method in Vorticity Formulation Before all, we show how the penalization method can be used successfully to model the flow of an incompressible fluid around an obstacle [START_REF] Angot | A Penalization Method to Take Into Account Obstacles in Incompressible Viscous Flows[END_REF]. 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. In this work, the penalization term is expressed using vorticity formulation. The main interest of the penalized vorticity formulation is that it replaces the usual vorticity creation algorithm in order to satisfy the no-slip boundary condition for vortex methods. This new technique avoids the convergence difficulties due to the particle creation on the solid boundaries (see [START_REF] Coquerelle | A Vortex Level Set Method for the Two-Way Coupling of an Incompressible Fluid with Collinding Rigid Bodies[END_REF] and [START_REF] Cottet | A Level-Set Formulation of Immersed Boundary Methods for Fluid-Structure Interaction Problems[END_REF]). We solve simultaneously the Brinkman equations in the solid and the Navier-Stokes equations in the fluid, considering whole the domain as a porous medium with zero (solid) or infinite permeabilities (fluid). Thus, defining the Reynolds number as Re = u re f l re f /ν, the non-dimensional penalized vorticity equation (or Brinkman-Navier-Stokes equation) reads

∂ω ∂t + (u.∇)ω = 1 Re ∆ω + ∇ × [λ χ S (u s -u)], (5)
where χ S denotes the characteristic function that yields 0 in the fluid and 1 in the solid, u s indicates the rigid body velocity which is zero in all this work since the body is fixed and λ = µΦ l re f /ρk u re f is the non-dimensional penalization parameter, in inverse proportion to the permeability of the medium (with k the intrinsic permeability, µ the viscosity, Φ the porosity of the porous material, l re f the height of the obstacle, ρ the density and u re f the reference 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. In addition it allows to compute the pressure as a continuous field on the whole domain including the solids. The zone variation is realized changing the penalization coefficient that defines the permeability of each region. Numerically, the fluid is considered as a porous medium with a very high permeability (λ = 0) and the bodies are considered as porous media with a very small permeability (λ = 10 8 ).

To discretize the penalized vorticity equation (Eqn. 5) in a vortex method, the equation is split in substeps. At each time step, one successively solves the unsteady penalization-convection-diffusion equations:

∂ω ∂t = ∇× (λ χ S (u s -u)) (6)
∂ω ∂t + (u • ∇) ω = 0 (7) ∂ω ∂t = 1 Re ∆ω. (8)
To solve Eq. 6 we use an implicit scheme ( [START_REF] Coquerelle | A Vortex Level Set Method for the Two-Way Coupling of an Incompressible Fluid with Collinding Rigid Bodies[END_REF]) and we set:

ω n+1 = ∇× u n + λ∆tχ S u s n 1 + λ∆tχ S . ( 9
)
where ∆t is the time step. The right hand side above is evaluated by second order centered finite differences. At this stage, grid vorticity above a certain cut-off is used to create particles in grid point surrounding cells and Eq. 7 is solved by a classical vortex-in-cell method [START_REF] Cottet | Vortex Methods -Theory and Practice[END_REF]. Particles are pushed with a RK4 time-stepping and are then remeshed on the original grid using the third order interpolation kernel M 4 [START_REF] Monaghan | Extrapolating b-Splines for Interpolation[END_REF]. Then, diffusion (Eq. 8) and Poisson equation (Eq. 4) are solved on the grid using FFT-based evaluations as described in [START_REF] Chatelain | A Fourier-Based Elliptic Solver For Vortical Flows With Periodic and Unbounded Directions[END_REF]. Grid values for vorticity and velocity are now available for time t n+1 and a new cycle of iterations can start. Moreover, the no-slip boundary conditions are naturally satisfied penalizing the vorticity transport equations.

This hybrid vortex penalization method has already been successfully used to simulate high Reynolds flow regimes past bluff bodies [START_REF] Cottet | A Vortex Penalization Method for Flows with Moving Immersed Obstacles[END_REF][START_REF] Cottet | A Vortex Immersed Boundary Method for Bluff Body Flows[END_REF].

Vortex Penalization Method for Solid-Porous-Fluid Media

This section is devoted to the study of flows in solid-porous-fluid media and aims to highlight the efficiency of the vortex penalization method to solve such problems. For three different solid, porous and fluid media configurations, it is of great importance to clearly understand the flow behavior at the fluid-porous interface in order to model the physics correctly. As described in [START_REF] Bruneau | Numerical Modelling and Passive Flow Control Using Porous Media[END_REF], we can consider five different flow regions from the solid to the free flow in the fluid. The first one is the boundary layer inside the porous medium, close to the solid wall. This one is very thin compared to the second region which is characterized by the homogeneous porous flow with Darcy velocity (numbers 1 and 2 in Fig. 1). In the vicinity of the porous-fluid interface, two transient layers can be recognized (numbers 3 and 4 in Fig. 1). The first one corresponds to the increase of the porous layer velocity reaching u i value at the interface and the second one to the fluid boundary layer standing from the interface to the free flow. The fluid boundary layer growth is then determined by u 0 -u i where u 0 denotes the velocity of the main fluid flow (number 5 in Fig. 1). The aim is thus to find out a way to solve the flow both 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 like avoiding to solve the porous flow enforcing appropriate porous-fluid boundary conditions [START_REF] Beavers | Boundary Conditons at a Naturally Permeable Wall[END_REF] or solving the governing equations of each region coupling Darcy equations and Navier-Stokes equations with a right treatment 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]. Nevertheless, these two approaches present some drawbacks since the first one neglects the porous medium physics and does not permit to have overall view of the problem and the second one is particularly difficult to handle especially because of the interface problem. The method presented here, based on the vortex-penalization technique, appears as a suitable alternative since it involves a unique equation (Eq. 5) for the whole domain and ensures an accurate modeling of each of the different regions thanks to the dimensionless penalization factor λ. Here, λ is expressed as λ = µΦ l re f /ρk u re f where l re f , ρ, u re f = 1 in this study and the porosity Φ is close to 1 as imposed by Brinkman equations [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 directly defines the different media according to the following equation obtained using implicit Euler scheme for the penalization velocity discretization:

u n+1 = u n 1 + λ∆tχ S . (10)
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. According to Eq. 10, the velocity in the fluid is not penalized ( u n+1 = u n ) and since λ = 0, the penalization term vanishes in Eq. 5, and we naturally recover the Vorticity Transport Equation (Eq. 3). On the contrary, the solid has a permeability coefficient k which goes to zero, it can be consequently modeled setting the penalization parameter λ to a very high value. In this study λ equals 10 8 in the solid, which vanishes the flow velocity in this region according to Eq. 10 ( u n+1 → 0) .

It has been proved in [START_REF] Angot | A Penalization Method to Take Into Account Obstacles in Incompressible Viscous Flows[END_REF] that solving Eq. 5 with such a value of λ was equivalent to solve Darcy's law in the solid. Furthermore, setting the λ parameter to an intermediate value, reasonably chosen between these two extreme values (λ = 0 and λ = 10 8 ), would model a porous medium in which the flow has a Darcy velocity u D (Fig. 1). For example λ = 1000 corresponds to a very low permeable medium and λ = 1 represents a very high permeability. The variation of λ corresponds to the variation of k that specifies the intrinsic porous material permeability. The accuracy and efficiency of the penalization method come from its capability to take into account these variations of λ and to capture the induced steep velocity variations at the different interfaces with a minimum number of discretization points.

Validation and grid convergence study

In order to verify the accuracy of the method, we present in this section a validation study of the computational method employed in this work. First we consider flows past semi-circular and circular cylinders at different regimes, and we compare the results with the ones of the literature. Secondly we present a grid convergence study for flow past a semi-circular cylinder, the obstacle of interest in this work. Each simulation presented in the following is performed considering the same geometrical setup: the obstacle under study has a dimensionless diameter of d = 1 and its center (in the case of a circular cylinder) or its back wall (in the case of a semi-circular cylinder) is located at (x, y) = (0, 0) in the computational domain D. The whole domain is meshed by an equispaced Cartesian orthogonal grid. As we use FFT-based evaluations to solve diffusion and Poisson equations, periodic boundary conditions are considered on the box walls and a correction of velocity is performed at each time step in order to satisfy the free stream velocity u ∞ = (u re f , 0) = (1, 0) imposed at the inlet.

Validation for flow past a circular and a semi-circular cylinder

The first part of this validation study deals with incompressible flow past a semi-circular cylinder.We compare our results to numerical and experimental results 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 length measurement, the Strouhal number S t (given by the shedding frequency) and the drag and lift coefficients, expressed by:

C D = 2F D ρu 2 re f d , C L = 2F L ρu 2 re f d . ( 11
)
where the drag and lift forces, respectively noted F D and F L , are computed according to the momentum equation given by Noca et al. 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]. The Reynolds number ranges between 100 and 300 in this study.

As can be seen in Table 1, 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 Fig. 2. Concerning in particular drag and lift coefficients, one can see from the left picture in this figure 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]. Indeed we find c D = 1.85 ± 0.15, cL = 0 ± 0.53 and S t = 0.230 to be compared with c D = 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]. The right hand side picture of Fig. 2 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 centerline in the x-direction where the streamwise component of the velocity vanishes. Again, our results agree rather well with the literature.

We complete this validation study by investigating the circular cylinder benchmark for which references addressing higher Reynolds numbers can be found in the literature. The Reynolds numbers considered in this part are the one we selected to perform the passive control study, that is to say Re = 550 (transitional regime) and Re = 3000 (highly transitional regime). The early-time evolution of the drag coefficient for both regime is reported in Fig. 3 c D ) and Strouhal number (S t ) for flow past a semi-circular cylinder at Re=100 and Re=300.

Grid convergence

We perform herein a grid convergence study of the vortex penalization method for flow past a solid semi-circular cylinder at Re = 550 and Re = 3000. This grid convergence is carried out on three grid levels. For Re = 550 the three consecutive mesh sizes are: h = 0.01, h = 0.005 and h = 0.0025 in the computational domain D = [-4, 8] × [-5, 5]. For Re = 3000, as the boundary layer is thinner, we perform the convergence study with three finer mesh sizes, respectively h = 0.005, h = 0.0025 and h = 0.00125. 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 results reported in Table 2, give the mean values of drag and enstrophy. For Re = 550, the results on the two finest consecutive grids are very close (Table 2). We can consider that the grid convergence is achieved and the adopted grid size for the further flow control simulations is h = 0.005 with ∆t = 0.0025. For Re = 3000 the grid convergence is achieved with h = 0.0025 and h = 0.00125. The corresponding studies will thus be performed using h = 0.0025 as the grid size and ∆t = 0.002 as the time step. 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. Table 2. Grid convergence study for flow past a semi-circular cylinder at Re=550 and Re=3000 through the mean values of drag coefficient ( c D ) and enstrophy ( Z).

Passive control

In this work, the vortex penalization method is applied to passive flow control around a two-dimensional semi-circular cylinder. The control is achieved covering this obstacle with a porous coating. The expected effect of a porous medium at the solid-fluid interface is to modify the vorticity generation of the boundary layer and the vortex shedding, and therefore to reduce drag forces and vortex induced vibrations improving the aerodynamic properties of the obstacle.

The subsequent flow control simulations are performed at transitional (Re = 550) and highly transitional regime (Re = 3000). As at these two regimes the flow is not turbulent, Direct Numerical Simulations (DNS) are performed to numerically solve Eq. 5. The semi-circular cylinder has a total dimensionless diameter of d = 1 including a porous layer of thickness τ, whose back wall is centered at (x, y) = (0, 0) in the computational domain D = [-4, 8] × [-5, 5] depicted in Fig. 4.

In order to analyze the effects of the control approach we will compare velocity, vorticity and pressure fields as well as global flow quantities like the drag force (F x ) and the enstrophy (Z) which measures the dissipation effects in the flow and defined by Z = D |ω| 2 dx. Note that in this study we will consider the drag force instead of the drag coefficient since the computation formula of the latter involves the diameter d of the obstacle which includes the porous layer.

Parametric study on penalization parameter λ for homogeneous porous layer

This section relies on the influence of the added porous layer permeability on the flow control behaviour and the efficiency of such a passive control. We perform a parametric study at Re = 550 and Re = 3000 considering four consecutive values of λ inside the porous coating, namely λ = 1 (high permeability), 10, 10 2 , 10 3 (low permeability) and compare the results to the fluid case (λ = 0) and the uncontrolled case. The fluid case means a semi-cylinder including a porous layer with permeability λ = 0, which corresponds to a smaller semi-cylinder with diameter d -2τ. In this section, the thickness of the layer is set to τ = 10%d = 0.1 for each case (Fig. 5).

As can be seen in Fig. 6, which represents dimensionless time history of global flow quantities at Re = 550, setting λ = 1 inside the layer clearly appears as the best solution in terms of flow regularization. Indeed, the mean value of drag At Re = 3000 the best solution for global flow regularization is also achieved setting λ = 1 inside the porous coating. This configuration leads to a drag and enstrophy reduction of respectively 6% and 30% compared to the uncontrolled case and one can see on Fig. 7 the important decrease of the signals amplitude. In terms of control effects, the main difference one can notice here in comparison to the parameter study performed at Re = 550 concerns the flow behaviour observed when λ = 10 inside the layer. Contrary to the other λ values, the flow field obtained with the latter is non-periodic (Fig. 7). This flow irregularity impacts the mean drag value, which is increased of about 15% compared to the uncontrolled case. We can summarize that except for a very permeable porous layer (λ = 1), a porous coating covering the whole face of the semi-cylinder is not a convenient tool to reduce the drag. It is only useful to regularize the wake. The problem seems due to the separation point that is located inside the porous medium. 6.2 Influence of the porous layer configuration on flow control efficiency This section aims to determine control devices which are both efficient in terms of flow regularization and easy to set up in terms of manufacturing constraints. According to the parametric study presented in the previous section, setting the λ function to 1 in the surrounding sheath enables important flow reduction effects. Nevertheless, this configuration turns out to be industrially difficult to implement because of the very high permeability of the added layer. On the contrary, intermediate permeabilities obtained fixing λ = 10 inside the porous interface are easier to handle but show quite poor improvements in terms of flow control, as reported in the parametric study. Taking into account these observations, we propose to compare the six test cases depicted in Fig. 8. Case 0 refers to the uncontrolled case, cases 1 and 2 correspond to homogeneous layer cases where λ respectively equals 1 and 10, and cases 3 to 5 refer to heterogeneous cases where the porous layer is split in two regions with different permeabilities. The simulations were performed taking the same numerical parameters as the one used for the parametric study. In a first place, the thickness of the layer is kept to 10%d = 0.1. The asymptotic mean values of global flow quantities and the time averaged vorticity, velocity and pressure fields are respectively given for each case in Table 3, and Figures 9 to 14.

First of all, the results reported here give further information concerning cases 1 and 2. They confirm that case 2 is not competitive as it increases the drag for Re = 550 and Re = 3000 compared to uncontrolled case (Table 3). They also highlight the efficiency of case 1 which induces a significant reduction of drag force (-27% and -6%) and enstrophy (-38% and -29%). These quantitative results are confirmed by the mean vorticity and velocity fields showing the smoothing of wake dynamics generated by the presence of the highly permeable layer. Indeed, the near wake structures are smaller and the back recirculation zone is sharply reduced (Figs. 9 and 12), implying an increase of downstream pressure (Figs. 11 and14) and thus a reduction of drag forces. Finally, the vortices swirl with lower velocity and for Re = 3000 the mean wake becomes axisymmetric (Fig 10 and13).

Let us now focus on the heterogeneous devices. For both regime, the case 3 shows benefits which are very comparable to those of case 1. We note that the only difference with case 2 is the presence of highly permeable poles in the layer. This change allows an eddy detachment from the wall with larger shedding structures, diminishing the back wall recirculation zone (Figs. 9 and 12) and a drastic increase of downstream near wall pressure for the Re = 550 regime (Fig. 11). The case 3 also enables to significantly regularize the flow in the wake with lower velocities, using highly permeable material introduced in both edges of the body (Figs. 10 and13). Nevertheless, the mix of high and intermediate permeabilities induces difficulties to industrially build such a coating and has no practical interest.

The cases 4 and 5 are more interesting for industrial applications. For these cases, the front part of the coating is no more permeable, but completely solid and porous layers are only placed in the two edge parts of the obstacle, which makes the practical implementation much more easier. Here, for both Reynolds numbers 550 and 3000, the enstrophy and the drag force undergo a tremendous decrease except for the case 5 at Re = 3000 where the drag force drop is limited to 7% because of the lower permeability of the porous layer (Table 3). Figures 9 to 14 confirm these results. They show a significant reduction of the mean velocity of the flow in the near and far wake (Figs. 10 and13) as well as the decrease of the crosswise dimension of the two counter-rotating vortical structures at the back of the obstacle (Figs. 9 and12). Also, the near wall pressure values (Figs. 11 and14), especially for the case 4, are sensibly diminished which is in agreement with drag results. These results emphasize the importance of the positions of the permeable zones that should be chosen considering the flow separation phenomenon as it preferably should not occur inside the porous material. Finally, test case 5, made of intermediate permeabilities on both edges of the body, corresponds to a more realistic device compared to the case 4 with very permeable edges. It can be considered as a suitable and affordable device for flow control.

Influence of the porous layer thickness τ on flow control efficiency

In the previous sections, the control effects were studied for different configurations using a constant porous layer thickness τ = 10%d. In this section, the effect of the thickness for different configurations will be investigated, and a special

Fig. 1 .

 1 Fig. 1. Velocity profile in the vicinity of a porous medium.

Fig. 2 .

 2 Fig. 2. (Left) Drag and lift evolution for flow past a semi-circular cylinder at Re=200. (Right) Normalized near-wake length L/d of a semicircular 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] (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).

Fig. 3 .

 3 Fig. 3. Drag evolution for flow past a circular cylinder at Re=550 and Re=3000. Comparisons of present method (solid red line) with Koumoutsakos & Leonard [16] (dashed line) and Ploumhans & Winckelmans [32] (solid black line).

Fig. 4 .

 4 Fig. 4. Computational domain used for passive flow control.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. (left) Uncontrolled case, (center) porous case, (right) fluid case.

Fig. 8 .

 8 Fig. 8. Cases 1 to 5 corresponding to different porous layer configurations. Case 0 depicts the uncontrolled case.

Table 5 .Fig. 15 .

 515 Fig. 15. Mean pressure profiles at the rear end of the body with a) τ = 10%d, b) τ = 5%d and c) τ = 2.5%d at Re=550.

Fig. 16 .

 16 Fig. 16. Mean profiles of the streamwise component of the velocity u at x = -0.05 with a) τ = 10%d, b) τ = 5%d and c) τ = 2.5%d at Re=550. d) Pattern of profile configuration at x = -0.05 showing the location of the profiles enlargement inside the layer proposed in a), b) and c).

c a s e 2 c a s e 5 cFig. 17 .

 2517 Fig. 17. Mean pressure profiles at the rear end of the body with a) τ = 10%d, b) τ = 5%d and c) τ = 2.5%d at Re=3000.

Table 1 .

 1 and one can see that the results are in excellent agreement with the considered references. Comparison of mean drag coefficient (

			Re = 100	Re = 300
	Authors	c D	S t	c D	S t
	Farhadi et al. [29]	1.73 0.180	1.955 0.22
	Present method	1.70 0.192	1.90	0.24

focus will be made for thin layers in order to reduce the use of porous materials and to find out a compromise between practical constraints and control efficiency. The same test cases (Fig. 5) are now considered setting the thickness of the layer to τ = 5%d = 0.05 and to τ = 2.5%d = 0.025. Table 4 and 5 respectively report the control performances in terms of mean global flow quantities, obtained for each case with τ = 5%d and τ = 2.5%d. For both regimes Re = 550 and Re = 3000, when comparing Tables 3, 4 and 5, a slight decrease of control effects is observed. This can be explained looking at the mean profiles of the streamwise component of the velocity, u at x = -0.05 plotted in Figures 16 and18. These profiles show higher acceleration of the flow field inside the layer when its thickness drops, inducing higher velocity gradients in the vicinity of the solid-porous-fluid interface and make the vortex shedding occurring earlier, resulting in a little reduction of beneficial effects. Nevertheless, the decrease of drag and enstrophy remains important for thinner porous layers, which enables to consider them as interesting passive control devices. Also, as shown in Figure 15, a direct correlation between the drag reduction and the back wall pressure forces is noticeable at Re = 550: when drag drops the back wall pressure rises. This fact is more unstable for the Re = 3000 as the flow is highly transitional (figure 17). Furthermore, for the later flow regime, decreasing the layer thickness in cases 1, 2 and 3 deteriorate the drag performances as the flow velocity slows down inside the thin layer around the body and the skin friction effects increase the drag forces. Now let focus on cases 4 and 5. One can observe in Tables 3,4, 5 that case 4 remains the best control device in terms of drag reduction, regardless of the coating thickness and Reynolds number. As discussed in the previous section, this is due to the fact that the front part of the obstacle is completely solid and that highly permeable zones are precisely located in the vicinity of the boundary layer separation point where the shedding starts. On one hand this configuration enables to efficiently decrease the shedding velocity of the flow damped by porous edges (see Figures 16 and18), on the other hand it prevents from the negative effect of porous coatings implemented in the front part of the semi-cylinder. For both cases 4 and 5, Figs. 16 and18 show for all thicknesses a uniform and relatively low acceleration of the flow field inside the layer. As the case 5 is the most practical benchmark (more than all other cases as the porous materials with moderate permeability and edge locations are easier to handle for the manufacturing of a real device), its control efficiency is important to be studied for all thicknesses: the control efficiency of case 5 is inversely proportional to the porous layer thickness τ. Indeed, in Tables 3, 4 and 5 one can clearly see that with τ = 2.5%d the drag reduction achieved with case 5, leads to -15% for Re = 550 and -10.5% for Re = 3000. This promising assessment seems to be due to the fact that the flow resistance inside a thin porous layer is less than a thick layer but the shedding is still well tuned thanks to the position of porous coatings in the body corners.

Re

CONCLUSION

In this work a passive control technique with porous coatings was implemented around a semi-circular cylinder in order to regularize flow and reduce drag forces for laminar and transitional flows. The direct numerical simulations were performed using vortex-penalization method which allows to compute simultaneously the three fluid-porous-solid media. In the numerical study, the effect of the porous layer permeability on the control results was first analyzed. Highly permeable media were more efficient, as expected. Then, several configurations of porous devices were compared. The results showed that introducing porous layers only at the top and bottom of the solid body permits to better manipulate the vortex shedding and to achieve relevant control performances. A further study on the influence of the porous layer thickness on the control efficiency was performed. It showed that a thin layer with intermediate permeability, introduced in both edges of the back wall, generates a significant drag reduction. As the ultimate target of this work is the drag reduction related to side mirrors, this type of configuration would be easy to implement industrially and at low cost, since it only consists in adding an intermediate-porous ring on the edge wall of the mirror. Although the present work needs to be complemented by threedimensional studies to confirm the quantitative control trends that were observed in the 2D experiments, it thus opens interesting perspectives towards efficient drag control devices.


Michel Jean

Jean Jacques

DYNAMICS OF ELASTIC OR RIGID BODIES WITH FRICTIONAL CONTACT: NUMERICAL METHODS

The mutual impenetrability of the parts of a mechanical system and their possible confinement by given external obstacles are unilateral constraints, the description of which results in joint inequalities arid equalities. The system dynarrùcs is thus governed by a set of very nonsmooth relations. The effect of dry friction at possible contacts mak:es evolution problems still less regular; velocity jumps may occur in this case, even in the absence of collision. In developing time-discretization procedures to solve these evolution problems numerically, the authors have chosen to face nonsmoothness without resorting to any regularization technique.

Introduction.

Let a mechanical system involve several parts, referred to as bodies. This lecture is devoted to computing the motion, when the possible contacts between some of these bodies, or between some of them and external obstacles with given motion, are subject to discussion at every instant.

Friction, affecting contacts, is a very intricate phenomenon. By restricting ourselves to the "dry" sort, we exclude the consideration of lubricants. Even so, the tribological data needed for an accurate description of contact interaction, when attempts are made at investigating engineering situations, are practically a:lways missing. Recall that, in usual instances, the alteration of contacting surfaces through wear and the accumulation of debris in the affected zone play an important role. Sorne authors also insist on the pressure-dependent crushing of micro-asperities and

introduce a specific nonlinear local compliance in the description of contact. An extensive discussion of dry friction, with a vast list of references, may be found in [l]. By representing friction through the traditional law of Coulomb, or a generalization of it to anisotropie situations, we choose to neglect these effects. Anyway, in our computing techniques of time-discretization, it is easy to make the parameters of the law of Coulomb depend, for each time-step, on some "status" variable determined by the past motion. More simply, at the beginning of each step, one may have the friction data depend, possibly in discontinuous way, on the available values of the system position and velocity.

Coulomb friction and unilaterality, i.e. the possibility of contact breaking at unknown instants, raise problems whose mathematical difficulty has long been recognized [START_REF] Lecornu | Sur la loi de Coulomb[END_REF]- [START_REF] Beghin | Sur l'indétermination de certains problèmes de frottement[END_REF]. Collisions, i.e. the sudden occurence of new contacts, causes still rougher irregularity.

A very accurate treatment of collisions would require data usually impossible to identify and to collect experimentally. The suddenly introduced contacts may persist in subsequent motion, or on the contrary, the colliding bodies may bounce against each other. The empirical coefficient of restitution bears little consistency beyond the special case of the collision of two otherwise unconstrained perfectly rigid bodies. An alternative concept has been proposed in [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF], avoiding in particular the energetical paradoxes pointed out by some authors [START_REF] Keller | Impact with friction[END_REF] [START_REF] Brach | Rigid body collisions[END_REF]; it is presented in Sec.6 below.

The traditional approach to the above problems consists of performing, at the first stage, some tentative calculation of the motion under the maximal regularity assumptions: persistence of existing contacts, no collisions, no occurence of singular points in the differential equations. From this calculation, one tries to identify the instants at which some contradiction arises: contact forces with inadmissible direction, collision, dynamical locking due to friction, ... At the first encountered of such signais, a special analysis of the instant behaviour (possibly with the help of infinitesimal time rescaling) is expected to yield the initial conditions for a possible subsequent phase of regular motion. Even in the no-friction case, this approach has long been recognized as requiring caution: the contacts which get loose after some of the critical instants are not necessarily those for which the previously calculated reaction exhibits inadmissible sign [START_REF] Delassus | Sur les liaisons unilatérales[END_REF] [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF].

The mathematical tools allowing for a global formulation of the concerned evolution problems, not based on the division of time into intervals of smooth motion, have been made available only in recent years. The existential study [START_REF] Jean | A system of rigid bodies with dry friction[END_REF] [15] and the corresponding numerical procedures [START_REF] Moreau | Numerical study of an oscillator submitted to unilateral constraints and dry friction[END_REF] are still under development.

For deformable bodies, the situation is still more complicated, since the contact status may, at every instant, vary from one region to another, delimited by unknown moving boundaries. So far, papers devoted to the computation of such dynamical problems are not very numerous [START_REF] Hughes | A finite element method for a class of contact-impact problems[END_REF]- [START_REF] Stein | Models of friction, finite-element implementation and application to large deformation impact-contact problems[END_REF]. In the event of a collision, shock waves are expected to propagate through the system; these are numerically considered in [START_REF] Chen | A new finite element technique for dynamic contact problems with friction[END_REF]. In the applications we have in view, the time taken by waves to travel the whole system is short enough for these to be inapparent in computation.

Anyway, it should be kept in mind that the presented models can only carry a crude representation of physical reality. This is a fact of life in all the applications of Science, but in the present topics, the lack of physical information is so patent that, however sophisticated the computation methods might be, seeking high precision in their implementation could be regarded as a waste of effort.

Analytical setting.

In the first Sections of this lecture, a mechanical system with finite number of degrees of freedom will be considered. Naturally, the numerical treatment of continuous media also relies on fmite sets of variables and we shall corne to this case in Sec. 9.

Let the possible configurations of the system be parametrized, at least locally, through generalized coordinates, say q = (q 1, q 2, ... , q n ). Such a reduction to finite freedom is assumed here to result from (bilateral) ideal constraints, namely the rigorous internai rigidity of the various parts and the possible operation of frictionless linkages.

After constructing this pararnetrization, one additionally takes into account some unilateral constraints whose geometric effect is expressed by a finite set of inequalities (2.1) fc/q) � 0, ae ( 1, 2, ... , K}, where f 1 , f 2 , ... ,f ic are given fonctions of q. Usually, such inequalities describe the mutual impenetrability of some parts of the system or the confinement of some of them by given extemal boundaries. Assuming, as we are doing here for simplicity, that the fonctions fa are constant in time involves that the possible extemal boundaries are fixed relative to the reference frame in use.

For brevity, parametrization is supposed scleronomic, i.e. time-independent. Then, for every imagined motion t�q(t) and for t such that the derivative q(t)e R n exists, the kinetic energy of the systems lets itself be expressed as a quadratic form in q , say

• 1 • i • j • (2.2)
'E/ q, q) = 2 A u (q) q q, w 1 th A u (q) = A / q). Then, as far as smooth , i.e. twice differentiable, motions are concerned, the system Dynamics is governed by Lagrange equations, written as an equality in R n (2.3) A(q) q = F(t, q, q) + r.

Expression F here comprises the terms of first differential order in Lagrange equations and the generalized components, relative to the parametrization (q), of some applied forces supposed given as fonctions of time, position and velocity. The element r = (r 1 , r 2 , ... , r n ) of R n is made of the generalized components of the totality of the contact forces corresponding to unilateral constraints.

Connecting r with the proper contact forces, which are vectors of physical space, resort to the standard machinery of Analytical Dynamics. Suppose that inequality f a. (q) � 0 expresses the mutual impenetrability of some pair of bodies '13 and '13'. Let us consider a position q 1 of the system such that equality f a. (q 1 ) = 0 holds, which means that in position q 1 the two bodies touch each other at some point of space called M a. . This we shall assume to be an isolated contact point, but other contacts, corresponding to different values of Cl, may also be effective between the same bodies in the considered position. Let t!denote the contact force that body '13 experiences at this contact point from body '13'; then '13' experiences from '13 the force -t!-. Defining the covariant components of this pair of forces, relative to the parametrization (q), rests on the kinematics of the respective particles of the bodies to which they are applied. One imagines a motion (usually qualified as virtual, nonnecessarily complying with the constraint inequalities) _ t�q(t) such that q(t 1 ) = q 1 . To every value, say ue R n , that the derivative q(t 1 ) may corne to take, correspond some values � and o/ � of the respective velocity vectors of the particles of '13 and '13' passing at point M a. . The corresponding power (usually said virtual) of the pair of contact forces then equals In ail the sequel, Greek indices, when repeated in an expression, shall not be understood as implying summation.

'li a. and 'Vd above are velocity vectors relative to the reference frame in use. The difference 'll a = 'V a -'Vd equals the relative velocity of body '13 with regard to body '13' at point M a. . In computation, this vector will be represented by its components relative to some chosen orthonormal base (n a. ,t a ,t � ), depending on q .We shall assume in the sequel that the unit vector " a. is normal to the contacting bodies at point M a and directed toward '13. Through the geometric investigation of the system, one expresses 'U a. from u in the form (2.4) where Gci: R n � R 3 is a q-dependent linear mapping. Such at least are the findings under the assumption we have made of a scleronomic parametrization. In general, Gci would also depend on t and a complementary term, constant in u, should be added to the right-hand side of (2.4).

By definition, the generalized components of the pair of contact forces make the element r a. of R n such that u.r a. = P a for every u in R n . Equivalently (2.5) This expression of r a through the transpose of the linear mapping Gci is known to remain valid in the absence of the scleronomy assumption we have made here for simplicity. In fact, the definition of generalized components, in Analytical Dynamics, rests on virtual motions at fixed time, i.e. chains of positions, say 't�q('t), depending on a real variable 't in differentiable way, with subsequent calculations made at constant t Observe that, at the present stage, the considered values of u are not required to be compatible with the contact geometry. This will be discussed in Sec. 3.

What precedes implicitely concerned a pair of bodies '13 and '13' which both were parts of the investigated mechanical system. The conclusion remains formally the same in the case where only '13 belongs to the system, while '13' denotes an external obstacle with prescribed motion. Among the two contact forces, only the force '1(:-experienced by '13 has to be kept into account in this case and one makes 'l. l a simply stand for the velocity vector � .

Summing up, under the smoothness assumptions made, the Dynamics of the system is governed by the differential equation ( 2.3) where the element r = I. r a is involved in the problem through (2.4), (2.5), u = q, and through a system of contact laws

(2.6) law( a, q, 'Ua , 'l(_ a ) = true.
The latter describe what, in physical space, happens at every contact a. In conditions (2.4) t<_> (2.6), the index a has to take all values such that f a (q) = 0 , since 'Ua, '1(:, etc. only make sense in the case of effective contact. Actually, in existential studies, as well as in computation, it proves convenient to make a range through a larger subset of { 1, 2, ... , K} provided the formulation (2.6) is contrived in such a way that f a (q) < 0 => 'i(:= O. In the same context, the matrix Gci which, strictly speaking, was only defined for q lying in the hypersurface f a (q) = 0 of R n has to be extended (in some smooth arbitrary way) to every q in a neighbourhood of this hypersurface.

One is looking for motions verifying the K inequalities (2.1) for every t. Instead of explicitely adjoining these to (2.3)-(2.6), we prefer to rely on some adequate formulation of the contact laws (2.6) for securing them.

Complete contact law.

DEFINITION. A contact law, i.e. a relation of the form (2.6), is said complete if it involves the three following implications

(3.1) (3.2) f a (q) < 0 => 'i(:= O. f a (q) � 0 => "a• 'l l a � O, (3.3)
The meaning of (3.1) is simply that the reaction force 'J(:' vanishes when contact a. is not in effect. Let us comment on the importance of (3.2). Put

(3.4) { { 'Ue R 3 : n a. .'U �O} if f a. (q) � 0 �(q) = R 3 if f a. (q) < 0,
called the set of the right-admissible values for the relative velocity of the two concemed bodies at point M a. . As in previous Section, the definition of n a. is assumed extended, so as to keep this notation meaningful for f a. (q) > 0, i.e. for system positions which violate the impenetrability of the considered bodies.

The following is easily established [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF].

PROPOSIDON. Let I be a time-interval with origin to and let a motion q:I�R n be defined through a locally integrable velocity function u:I�R n by

(3.5)
t � q(t) = q(t o ) + r u(s) ds.

to If 'll a. (t)=G� u(t) belongs to �(q(t)) for almost every t and if inequality fi q( t) )�O holds at the initial instant t 0 , then this inequality holds for every te I.

In other words, provided the initial position is correct, the impenetrability condition f a. �O is automatically taken care of by (3.2). Observe that this proposition is sensitive to the ordering of time. In the symmetric statement, involving instead of the initial instant t 0 , some final instant, one should replace � by -� , which may be viewed as the set of the left-admissible values of 'U a. . The importance of (3.3) will only become apparent in further Sections, devoted to numerical algorithms and to the study of collisions.

From our standpoint, synthetic formulations of contact laws, possessing properties (3.1)-(3.3) without a priori statement of them, should be preferred. The theoretical or numerical developments based on such formulations offer a securizing consistency and should accept with better robustness the further adjustments needed to improve the agreement with physical reality. Such a synthetic formulation for the frictionless contact law has been presented in [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF]. Several ways of doing the same for dry frictional contact will be developed in further publications.

Coulomb contact.

In this lecture, the status of complete contact law will be given to Coulomb contact by asserting from start that properties (3.1) to (3.3) hold. Therefore � can be nonzero only if f a.

(q)�O and n a. . 'll a. = 0, and there just remains to develop this case.

In some early publications of one of the authors [START_REF] Moreau | Sur les lois de frottement, de plasticité et de viscosité[END_REF]- [START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF], the traditional statement of Coulomb law, and a plausible extension of it to possible anisotropie friction, were replaced by equivalent assertions involving some concepts now standard in Convex Analysis. These statements are listed below under the title of "First formulation". Their advantages were shown to lie in possible handling through minimization techniques and in the simplification brought to the treatment of mechanical systems combining several frictional contacts.

In what follows, the data for every contact point consist of a closed convex cone C, in the linear space of physical vectors, called the Coulomb cone of the considered contact. Let us use the base (n, t, t') introduced in Sec. 2, omitting for brevity the contact index a . Vectors t and t' generate the two-dimensional linear subspace T, called the tangent plane at the considered contact. The Coulomb cone is assumed to contain n in its interior and to lie in totality on the same side of T.

Giving the cone C is equivalent to giving its section at unit distance of T or, more conveniently, the orthogonal projection of this section onto T, say '1> 1 .This is a closed convex subset of T, containing the origin as interior point. Giving '1> 1 is in turn equivalent to giving its supportfunction. One preferably uses the image of this function under the symmetry relative to the origin, i.e. It will be seen in a moment that <1> 1 may be called the dissipationfunction at unit normal component.

For the traditional isotropie Coulomb friction, '1> 1 consists of the disk centered at the origin, with radius equal to the friction coefficient y, so <1> 1 ( 'U) = y 11 'U Il .

FIRST FORMULATION OF COULOMB LA W. One decomposes the contact force in the form It is first stipulated that the normal component p is nonnegative. As soon as this real number is treated as known, one may denote by <l> the sublinear function p<J> 1 , by '1J the set p'1> 1 and by 'l'p the indicator function & this set (namely 'I' =O in '1J P and P 'l' _p =+ 00 elsewhere). Then Coulomb law lets itself be expressed by �y of the equivalent following relations between the elements 'U and � of T: Recall that a'I' , the subdifferential of the indicator fonction of the set '1J is nothing but the (oufward) normal cone to this convex set at the considered poin�. In Convex Analysis, the normal cone to a convex set is defined at any point, say x, of the embedding linear space. This cone is empty if and only if x does not belong to the set ; otherwise, it contains at least the origine and may reduce to this (in particular when x lies in the interior of the set). So (4.3) involves that �E '1J P .

As a point of vocabulary, one may express (4.4) by saying that, for fixed p, the (multivalued) relation between the velocity 'l1 and the frictional force � adroits the convex real fonction <l> p as pseudo-potential.

The equivalence of (4.3) to (4.4) is due to the fact that the fonction 'I' (the indicator of '1J P ) and the fonction 'V� <1> (-'V) (the support fonction of '1J ) Jake a pair of FencJiel conjugate functions. 'rhen 1{. T and -'ll are conjuga�e points relative to this pair of conjugate fonctions, a situation classically expressed by Fenchel equality. Conditions (4.5) are merely the special form that Fenchel equality takes in the present context; the equality on the right justifies the name of dissipation function (for the value p of the normal component) given to <1> . Conditions (4.6) amount to rewriting (4.5) in the form of a "principle" of maximafdissipation.

In short, when p is treated as known, Coulomb law is exactly similar to a law of perfect plasticity, well known to generate mathematical problems in the form of variational inequalities. In situations involving frictional contact, it is commonly found that the determination of p is in turn governed by a variational inequality involving 1{. T as parameter, so one is in the presence of two coupled variational inequalities. The precise proof of this equivalence may be based on standard arguments of Convex Analysis [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF]. Let us only comment on (4.7) by observing that it exhibits the expected features. The subdifferential � ç (1Ô is nothing but the normal cone to the convex set C at point � Since this cone 1s empty for 1{.<l C, relation ( 4. 7) forces 1{. to belong to C whichever is 'l1 in T. If 1{. lies in the interior of C, the normal cone reduces to the set { 0}; hence it is only for 'U = 0 that 1{. can take values interior to C. The normal cone to C at the origin consists of the totality of the polar cone c 0 of C. In most usual situations '1J 1 is bounded; then through standard propositions of Convex Analysis, C 0 adroits -n as an interior point Consequently, the projection of C 0 onto T equals the whole of T ; in other words, ( 4. 7) lets the value 1{_ = O compatible with any 'l1 ET. o e <hv/'lô + ae('U).

SECOND FORMULATION OF COULOMB

The latter has been successfully used in the numerical treatment of elementary three-dimensional contact problems [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF]. Expressions based on standard Convex Analysis and involving, as (4.9) does, two independent "pseudo-potential" functions could also be applied to the formulation of other resistance laws, e.g. plasticity laws with flow rule "non associated" with the yield criterium. In the same line, a new type of numerical procedures of the variational sort has recently been proposed by G. De Saxcé and Z. Q. Feng [START_REF] De Saxcé | New inequation and functional for contact with friction: the implicit standard material approach[END_REF] under the name of the "method of bipotentials". This method may also lead to alternative statements of Coulomb law.

Collisions.

When the unilateral constraints of impenetrability are taken into account, the treatment of dynamical problems cannot in general exclude the possibility of collisions. Unfortunately, in practical instances, the physical data required for accurately predicting the consequent motion are missing. The assomption of elastic shocks is commonly made, which means that the considered collisions involve no energy loss. This assumption yields a scalar equation, which anyway is not enough to predict the after-shock velocities of all system parts, if several unilateral contacts are present at the collision instant.

Energy preservation can rarely be justified. The concept of a collision between bodies which, in the rest of the motion, are treated as perfectly rigid can only be the schematic representation of a "very short" interaction involving "very large" contact forces. Even if it is safely established that, in spite of the large stresses involved, the material of which the system is made behaves elastically, the collision locus is sure to be the source of dynamical disturbances which propagate through the whole system. If the system is linked to supports, the se disturbances will also affect the external world, carrying away an irrecoverable amount of energy. In the case of a bounce, various parts of the system will still be vibrating when the colliding bodies get loose from each other. By considering all parts as rigid in the subsequent evolution one treats vibration amplitudes as "microscopie", but this does not imply that the vibrational energy is negligible.

Of course, the possible inelastic deformation of the colliding bodies, however small its amplitude is estimated, may also cause energy loss.

What precedes holds true even if, during the "very short" duration of the collision, the contact of colliding bodies is considered as frictionless. Things become still worse if friction is ta.ken into account. Drawing the momentum balance (or performing the equivalent calculation in the setting of Analytical Dynamics) would require the knowledge of contact percussions, i.e. the integrals, over the "very short" duration of the collision, of "very large" contact forces. Even if one makes the hazardous assumption that, under such high contact pressures and rapidly varying conditions, friction obeys the same laws as in regular situations, there remains to observe that Coulomb law consists of a nonlinear relation between forces and velocities. Commuting it with time-integration, so as to obtain a relation between the contact percussion and the average contact velocity, is thus illegal.

A favourable special case should however be mentioned: if during the collision, the relative velocity 'l1 of the contacting bodies (see Sec. 4 above) has constant direction, i.e. 'll(t) = cr(t) 'l1 0 , where 'l1 0 is a constant vector and cr(t)>O, Coulomb law yields that the contact force 2«t) remains in a fixed closed convex cone, namely the set of values of 1?.. that (4.7) associates with 'l1 0 . Therefore, the contact percussion also belongs to this cone, so it appears to be associated with 'l laverage = cr average 'llo• An early attempt at analysing the variation of 'l1 during a frictional collision through rescaling in time, was made by G. Darboux [START_REF] Darboux | [END_REF] and later improved by J. Pérès [START_REF] Pérès | Mécanique Générale[END_REF]; see also [START_REF] Keller | Impact with friction[END_REF]. This consisted in deriving from the equations of Dynamics, joined with Coulomb law, a system of first order differential equations governing the fonction t � u(t) during the collision, while the position q was treated as a constant. Sorne assumption had to be added in order to decide the end of the process. In this approach, the bodies were considered as perfectly rigid. Physical situations likely exist, where the assumptions made are valid. The cruder treatment we propose in the sequel may be equally usefol in practice.

The dissipation index.

A collision, occuring at some instant t s , is expected to generate a discontinuity for the velocity fonction t � u(t). Traditionally, the left-side limit u-(t s ) and the right-side limit u + (t s ) are assumed to exist. A consistent mathematical framework where this existence is secured, is provided by assuming that u is a fonction of bounded variation [START_REF] Moreau | Bounded variation in time[END_REF]. Theoretical arguments in favour of this assumption are developed in [START_REF] Moreau | An expression of classical dynamics[END_REF].

Let us apply the usual trick of integrating both members of (2.3) over the "very short" duration of the shock. The position q(t) is treated as a constant and the term F remains bounded, so it yields a negligible integral. One thus obtains (6.1) where the element s of R" equals the time-integral of r. Force vectors of physical space are connected with their analytical representatives through relations (2.5). Like q, the matrices a:t are treated as constant during the shock, so one has (6.2) where .11-, the time-integral of !!( a , may be called the shock percussion at contact ex.. During the shock process, the relative velocity 'U a of the contacting bodies rapidly varies from � = Gci u-(ts) to � = Gci u + (ts). Now cornes the embarassing question of connecting the unknown vectors S a with the sought motion. Rather than attempting to calculate the evolution of 'U a during the shock, we propose, as a first approximation, to accept that a contact law similar to (2.6) holds between S a and a certain weighted mean m 1-ô _ l+ô

(6-3) 'U a =2 'U a + y � '
where 8 denotes a real constant, called the dissipation index [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF].

Sorne confidence in this procedure may be gained by drawing the energy balance of the collision. From the expression (2.2) of the kinetic energy, orte readily derives the jump of this function in the form

'E + -'r = -2 1 A .. (u� -u-:-)(u� + u-:-). C C IJ 1 1 J J
Here Ai j stands for the matrix A i /q(ts)), so (

E;-'Ç = f (u + + u-).s . 6.1) yields '
By making use of (6.1), (6.2) and ( 6.3), one obtains the signed decrement of kinetic energy

J J a a

As far as �O. the first term on the right-hand member is nonnegative ; we propose to call it the structural dissipation. The nonnegativity of 8 may be viewed as asserting that no external supply of energy cornes to enhance the possible bounces. The last term depends on the adopted contact law and may be called (with the minus sign included) the frictional dissipation.

The following shows the importance of including in contact laws the "completeness" assertion of Sec.3. In the present context, (3.3) becomes ��o � " a •'U: � o.

In view of the definition (6.3) of 'U:;, the latter inequality yields l+ô + o-1

(6.5) 2 " a • 'U a � 2 " a • 'U ; .
For an effective collision, "a•� is strictly negative. If ô was strictly greater than one, " a • 'U� would therefore be strictly negative, in contradiction with the impenetrability constraint.

In short, the parameter o must be assigned a value in the interval [0, 1].

Making now use of (3.2), one concludes that, provided s'1-'# 0, (6.5) actually holds as an equality, hence + 1 --<> n a: 'U a.

= -1 +o " a. • 'U a. , which allows one to identify the expression 1--<> /1 +o with the coefficient of restitution, traditionally introduced for the description of two-body collisions. But when studying cases where several contacts are in effect at the collision instant, it is essential to observe that this equality does not necessarily hold if s'1-= O. An elementary example of such a situation is provided by the rocking oscillations of a column in contact with a fixed horizontal plane.

Pri n ciple of algorithms.

Let [t 1 , tF] , tp=t 1 +h, denote one of the intervals of a time-discretization. Sorne approximate values q 1 , u 1 of the fonctions q and u at instant t 1 result from the preceding step. In order to calculate an approximate value of u at tF, say uF , one discretizes (2.3) into (7.1) where u L = u 1 + A-1 Fh , the "loose velocity", is the value one would find for uF in the absence of contact force. Here, A, F, a ;!* denote the values calculated for the corresponding fonctions at a point (t, q, u), chosen as typical of the time-step and computed in some explicit way. In contrast, the procedure we shall use to calculate uF may be qualified as implicit. In harmony whith the notations of the preceding Section, s'1denotes the integral of the contact force '.K!on [t 1 , tF]. Therefore the algorithm will be able to treat possible shocks on the same footing as smooth motions.

The decision of applying a computing policy of the implicit type, induces to invoke the contact law, for every selected a, as a relation of the form (2.6) holding between s'1and the final local velocity, namely in view of (2.4)

(7.2)
After applying a i to both members of (7.1) and (7.2), one obtains (7.3) 'Up F = 'Up L + f a i A -l a ;!* .s'1with 'lip L = a i u L , to be joined to the contact law (7.4) law( a, q, 'Ua. F , .s'1-) = true.

Of course a and � will be restricted to range in the (possibly empty) subset ACT of { 1, 2, ... , K} corresponding to the contacts considered as active in the rime-step. The simplest criterium for determining this set is based on the signs that the fonctions fa take at point q 1 ( or at some other test position). Solving (7.3)-(7.4) will yield the vectors .s°', then Up through (7.1). After that, the computation step may be concluded by taking 1 q F = q l + 2 h (u l + U p).

As soon as several contacts are active, (7.3)-(7.4) can only be treated by successive approximations. Let us stress that, in general, one cannot expect uniqueness for the solution vectors §1-. Dry friction commonly generates hyperstatic problems. In favourable cases, the calculation may however result in a unique u F The technique we have applied in ail recent experiments is a relaxation method, based on the solution of the single contact problem. In particular, for the two dimensional case, the case of a single contact can readily be solved, yielding for every a an easily progràmmedpiecewise affine mapping 'll a.1., �.11-. Starting with a tentative set of values for the vectors .s°', ae ACT, one makes a cyclically run through ACT. Each time, one updates the corresponding approximate value of .s°' bX solving a single contact problem, where the previously obtained approximates jl, �:;éet, are treated as known forces applied to the system. The above computation is valid even if a collision occurs in the rime interval [t 1 , tp], a situation revealed by the appearance in ACT of some contacts which were not active at the antecedent time-step. Clearly, by invoking 'U aF as value for the local velocity in (7.4), we make the procedure treat possible collisions as inelastic, i.e. the dissipation index equals one. By introducing into (7.4) the same weighted mean as in (6.3), one could impose a different value of 8. It proves more expedient to program the algorithm for the case 8 = 1 only, getting a result one may denote by u F = solver(u L ). Due to the fact that relations (7.2)-(7.4) are positively homogeneous with degree one, the time-step result for arbitrary 8, is found to be

u F = [2 solver(u L ) -(1-8) u L ]/(1 +8).
It has long been recognized [START_REF] Lecornu | Sur la loi de Coulomb[END_REF] that in dynamical problems involving Coulomb friction, some "catastrophes" may occur, generating velocity jumps in the absence of any collision. The above algorithms work all right in such situations [START_REF] Moreau | Unilateral Contact and Dry Friction in Finite Freedom Dynamics[END_REF].

In most our applications, the end of the iteration procedure is not decided upon the quasi-stationarity of the sequence of approximants, but on testing the precision at which the required relations are satisfied. In all experiments, the process has behaved properly, even for unrealistically large friction, but we have no proof of its convergence.

Deformable bodies.

The deformability of the considered bodies is accounted for through the introduction of a finite number of additional variables. In some instances these variables may arise from a modal representation. More frequently, the displacement fields, with respect to some reference position 2 are discretized by means of finite element meshes. Let the real variables q 1 , q , ... , q n refer to the corresponding degrees of freedom. Then the displacement at time t, say Q(t, x) e R 3 , of a particle x belonging to one of the bodies, admits the approximate representation Q(t, x) = q i (t) Ei(x), where E 1 , ... ,En denote R 3 -valued functions associated with finite elements. In the course of a motion t�q(t), the approximate velocity v�ctor of x is obtained by time derivation, as a linear expression in the derivatives q 1 • Let a n denote a material boundary part, that preliminary investigation has e detected as possibly entering into contact with another body of the system or with some external obstacle. For every particle xe a en, one may express the relative velocity of x, with regard to its potential antagonist, in the form Here 'U i and G� have been constructed in such a way that the value Ue R 3 consists of the c6mponents of the said velocity with respect to a x-dependent orthonormal base made of vectors normal and tangent to the antagonist body at the expected contact point. Index i has only to run through a subset that we shall suppose to be { 1, ... , p}, corresponding to the finite elements involved in the representation of aen.

In the principle, a contact law of the form (2.6) should be stipulated pointwise throughout a en, enforcing a relation between the local relative velocity and the density of reaction force R(t,x). Actually, one has to rely on an approximate representation of this force density, using an adequate set (F.) of R 3 -valued base functions, in the form . R(t,x) = 2S (t).F / x), with j runriing in ( l , ... ,p} and 1SE R 3 . One constructs the matrix P .. = fa ,... E.(x).F.(x) da (da= area), l j C U 1 j so the power of a distribution of reaction forces is expressed as (P Gq q ( t) ).1{{ t) or also (P 'U(t)).1{!:t), where 1( and 'U denote two collections of p elements of R 3 .

Contact is discussed by introducing, instead of the functions f a ' the normal gap qN(t, x), between particle x and the antagonist boundary. Then the contact laws adopted in the preceding Sections yield the classical Signorini condition which, for every XE a en, may be written as (8.2) RN(t, x)�O and 'v'sNe R + : 'JN (t, x)(sN-RN(t, x)) �-Now � (t, x) lets itself be represented in the form Q Ni (t) E i (x). Let us denote by P T and P N the blocks of the matrix P referring to tangential and normal components. By integrating both members of this inequality on êlé2, one obtains

The same technique may be applied to Coulomb law, invoked in the form (4.6). For every j in { l , ... ,p}, a convex subset 'IJ of the tangent plane is introduced as in Sec.4. The multiplication of these convex sets by the respective nonnegative normal reactions, yields a convex subset of the product space, denoted by � 'IJ. Tuen (4.6) becomes

The numerical treatment of the problem is thus quite similar to what has been presented in the foregoing for the case of a system of rigid bodies.

The base vector fonctions F. through which one represents the contact force densities are always chosen witfi nonnegative components. If the finite elements used in representing displacements are, for instance, of type T3 or Q4 the components of fonctions E. are nonnegative, as well as the elements of the matrix 1 P. In such a case, the approximation of the contact forces distributed on êl é2 by forces applied to the nodes proves to be a valid computation technique. On the contrary, the use of finite elements of type T6 or Q8 yields components of E i with no definite sign and then the above machinery provides much better results than the approximation by forces concentrated at the nodes [START_REF] Jean | Unilateral contact with dry friction: time and space discrete variables formulation[END_REF] [START_REF] El Youssoufi | Contact unilateral avec frottement sec en milieux continus discrétisés[END_REF].

In the event of a collision, an important differènce should be stressed between the treatment of perfectly rigid bodies and that of deformable ones. In the latter case, the total mass of the material which direct/y takes part in the contact is zero. Since this material has by itself no inertia, any condition distinguishing beween its right and left-velocities, like those stipulated in Sec. 6 for systems of rigid bodies, is dynamically meaningless. Signorini conditions (8.2), joined with friction laws and, of course, the dynamical equations satisfied throughout the deformable bodies, are all what govems collisions.

Examples.

The first figure presents computed successive positions for a two-dimensional system of 18 rectangular rigid bodies submitted to gravity. This is meant to model a pillar made of blocks superposed without mortar. As a crude simulation of earthquake, the horizontal ground is given a motion of sinusoidal horizontal translation. The same value 0.3 has been adopted for the Coulomb friction coefficient between blocks and with ground. The dissipation index equals 1.

When two blocks touch each other along a line segment, interaction is treated as concentrated at vertices (as if block edges were slightly concave). Even so, the problem of calculating contact forces at every instant is likely to be "hyperstatic" and the motion consequent to initial data is not expected to be uniquely determined. Calculation is also very sensitive to constants: if the height of some block is slightly different from that of its neighbours, the motion is visibly altered. However, all simulations agree with a consistent ruin pattern: the lower blocks progressive/y get apan, creating a gap.

The second figure shows successive positions of two elastic blocks, the deformation of which is accounted for through two-dimensional finite element meshes, as developed in Sec. 9 above. The two blocks are made of the same material (E = 1300 MPa, v = 0.2, density = 7.8 10 3 KgJm\ The larger block is fixed along its lower edge and initially at rest. Before hitting it, the upper block, with natural size of 0.04 mx 0.04 m, has a motion of translation with velocity components (-75 m/s, -50 m/s). Friction coefficient equals 0.1. Ten successive positions are shown, from t = 0 to t = 0.75 10-3 s. After a phase of sliding contact, the overall motion results in a bounce, with persistent vibrations of each body.

The two series of drawings are extracted from microcomputer animations which allow for a more vivid perception of the effect of contacts.

( 4 . 1 )

 41 'U-) <j> 1 ('U) = sup{ 'U.'TE R : 'Te-'1> 1 }.

( 4 . 2 )

 42 � = � + p n , with pe R and n. � = O.

  E a\j/p(�) -� E aq>p('U) � E '1Jp and -'U.� = <j>p('U) � e 'D P and V'Te '1J P : -'U.'T� -'U.� .

  LAW. A formulation avoiding the separation of the normal and tangential components of 1{. has more recently been proposed [7][26]. It involves the indicator fonction 'If c of the Coulomb cone, a convex lower semicontinuous fonction defined on the three-dimensional space of physical vectors, with values O or + oo. Assuming again that (3.1)-(3.3) hold, what precedes is found equivalent to (4.7) -'U E proj T a'lf / �).

  An alternative statement to (4.7) is obtained by defining on the three-dimensional space of physical vectors, the fonction 9 with values in [O, + 00 )

( 6 . 4 )

 64 'r -'E + = -2 1 ô A.. (u7' -u-:-)(u� -u-:-) -:E ef? . �-

(8. 1 )

 1 U(t, x) = 'LJ.(t) E.(x) = Gq q• j(t) E.(x).

V v.,., .... �;;: y;;;J7.,.,.,.,�;;: 9.,.,.,., V Horizontal ground oscillation: period 0.4 s., range 15 cm. Pillar initial size: width 90 cm, height 360 cm. Successive states shown at one period intervals.

Friction coefficient between blocks and with ground: 0.3

Il FRICTIONAL COLLISION OF TWO ELASTIC BLOCKS Larger block initially at rest and fixed along its lower edge.

Smallcr block collides wilh zero angular velocity and velocity componcnts -75 m/s, -50 m/s.