I. INTRODUCTION

In recent years, a significant research effort has been dedicated to the development of new flow energy-harvesting systems and the characterization of their efficiency. Many of such systems are based on flow-induced vibrations of structures, such as vortex-induced vibrations [START_REF] Taylor | The Energy Harvesting Eel: A Small Subsurface Ocean/River Power Generator[END_REF][START_REF] Bernitsas | VIVACE (Vortex Induced Vibration Aquatic Clean Energy): A New Concept in Generation of Clean and Renewable Energy From Fluid Flow[END_REF][START_REF] Grouthier | Energy harvesting by vortex-induced vibrations in slender structures[END_REF], flutter of wing profiles [START_REF] Erturk | On the energy harvesting potential of piezoaeroelastic systems[END_REF][START_REF] Bryant | Aeroelastic flutter energy harvester design: the sensivity of the driving instability parameters[END_REF][START_REF] Xiao | A review on flow energy harvesters based on flapping foils[END_REF], of cylinders in axial flows [START_REF] Singh | Energy harvesting from axial fluid-elastic instabilities of a cylinder[END_REF][START_REF] Singh | The effect of non-uniform damping on flutter in axial flow and energy harvesting strategies[END_REF], or flags [START_REF] Tang | Cantilevered flexible plates in axial flow: Energy transfer and the concept of flutter-mill[END_REF][START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Akcabay | Hydroelastic response and energy harvesting potential of flexible piezoelectric beams in viscous flow[END_REF]. The present work focuses on the latter, and considers a flag covered with piezoelectric patches that convert mechanical deformation into electric charge transfer inside an output circuit.

In this article the term flag refers in fact to a clamped-free plate in an axial flow. It is well-established that such structure can undergo strong self-sustained oscillations once a critical flow velocity is overpassed (see the recent review by Ref. [START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF] and references therein). These oscillations originate from an instability of the equilibrium position, resulting from the interaction of the flow forces with the solid's inertia and rigidity. This instability can be predicted from the linearized dynamics equations and different methods and models have been proposed depending on the aspect ratio of the flag, ranging from two-dimensional models (large span, [START_REF] Kornecki | On the aeroelastic instability of two-dimensional panels in uniform incompressible flow[END_REF][START_REF] Guo | Stability of Rectangular Plates With Free Side-Edges in Two-Dimensional Inviscid Channel Flow[END_REF]) to slender body models (short span, [START_REF] Lemaitre | Instability of a long ribbon hanging in axial air flow[END_REF]). The general case of a three dimensional fluid flow around a rectangular plate showed the limits of applicability of both models [START_REF] Eloy | Flutter of a rectangular plate[END_REF][START_REF] Doaré | Flutter of an elastic plate in a channel flow: Confinement and finite-size effects[END_REF]. Linear analysis is however not able to predict the amplitude of flapping in the saturated regime, hence it is necessary to consider the nonlinear dynamics of the flag to be able to determine how much energy can be harvested. Beyond full numerical simulations of the coupled fluid-solid problems, simplified fluid-solid nonlinear models have been proposed in the limit of two-dimensional flows [START_REF] Alben | Simulating the dynamics of flexible bodies and vortex sheets[END_REF][START_REF] Michelin | Vortex shedding model of a flapping flag[END_REF] or slender body problems [START_REF] Singh | The effect of non-uniform damping on flutter in axial flow and energy harvesting strategies[END_REF][START_REF] Eloy | The origin of hysteresis in the flag instability[END_REF].

Piezoelectric patches have traditionally been studied to couple structural mechanics to electrical circuits in the context of passive damping of structural vibrations [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF], of active control of vibrations [START_REF] Preumont | Vibration Control of Active Structures: An Introduction[END_REF] and of course energy harvesting [START_REF] Lefeuvre | A comparison between several vibration-powered piezoelectric generators for standalone systems[END_REF][START_REF] Anton | A review of power harvesting using piezoelectric materials (2003-2006[END_REF][START_REF] Elvin | Advances in energy harvesting methods[END_REF][START_REF] Chen | Study of a piezoelectric switching circuit for energy harvesting with bistable broadband technique by work-cycle analysis[END_REF]. The use of a piezoelectric flag to convert the kinetic energy of a flow was initially introduced by Ref. [START_REF] Taylor | The Energy Harvesting Eel: A Small Subsurface Ocean/River Power Generator[END_REF], to produce energy from the flapping of a flexible membrane forced by the wake of an upstream obstacle. In the last five years, several experimental, numerical and theoretical studies have been conducted to study the possibility to exploit the flutter instability itself using piezoelectric plates [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Dunnmon | Power extraction from aeroelastic limit cycle oscillations[END_REF][START_REF] Akcabay | Scaling the dynamic response and energy harvesting potential of piezoelectric beams[END_REF]. In these works, the plate was either entirely or partially covered by one piezoelectric pair [START_REF] Dunnmon | Power extraction from aeroelastic limit cycle oscillations[END_REF][START_REF] Akcabay | Scaling the dynamic response and energy harvesting potential of piezoelectric beams[END_REF], or completely covered by a large amount of small piezoelectric elements, so that a continuous model for a homogenous plate can be considered [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. In the latter case and for purely resistive output circuits, it was shown that the maximum efficiency is obtained for large mass ratio when the timescale of the electrical circuit is equal to the fluid-solid instability timescale [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. For practical applications, however, only a finite number of piezoelectric elements can realistically be used. The goal of the present article is therefore to provide some understanding on the impact of such a discrete coverage on the dynamics and performance of the system, a critical question that has remained so far unanswered.

In contrast with simplified representations of the energy harvesting process (e.g. pure damping, [START_REF] Tang | Cantilevered flexible plates in axial flow: Energy transfer and the concept of flutter-mill[END_REF]), such an approach provides a fully-coupled description of the fluid-solid-electric problem, taking into account both the transfer of energy from the solid motion to the output circuit and the feedback effect of the piezoelectric patches on the flag's dynamics. The nonlinear dynamics of the flag (even in the absence of any piezoelectric coupling) exhibit strong inhomogeneities of the local deformation and motion along the plate: near the clamped leading edge, the flag has little displacement but significant curvature, while near the free trailing edge, the flag has zero curvature but undergoes large displacements. Local stretching and compression of the piezoelectric patches resulting from the flag deformation are responsible for the charge transfer in the output circuit; the location along the structure of the device converting mechanical to electrical energy is therefore expected to strongly affect the efficiency of the system. A similar question was recently addressed in Refs. [START_REF] Singh | Energy harvesting from axial fluid-elastic instabilities of a cylinder[END_REF][START_REF] Singh | The effect of non-uniform damping on flutter in axial flow and energy harvesting strategies[END_REF] in the case of a pure damping model.

The present work focuses on the optimization of the distribution of a small set of piezoelectric elements (one, two or three pairs) positioned on a plate in an axial flow. In section II, the physical model and equations are presented for the fluid-solid-electric system with a finite number of discrete piezoelectric patches. A method of simulation based on the weakly nonlinear form of these equations is presented that is appropriate to obtain the nonlinear dynamics near the instability threshold. The influence of the number of piezoelectric elements on the stability of the system is investigated in section III. In section IV, the impact of this arrangement on the nonlinear dynamics and energy harvesting efficiency is studied using experiments and numerical simulations. Finally, a parametric numerical study is performed to determine the optimal arrangements in the case of one, two or three piezoelectric elements.

II. PHYSICAL MODEL

The system considered in this work consists of an elastic flexible plate of length L and width H immersed in an axial flow of uniform density ρ and velocity U ∞ . The plate is inextensible and clamped at its leading edge. For simplicity only purely two-dimensional motions of the plate are considered. The plate's geometry is therefore completely determined by the local orientation of the plate with respect to the flow, θ(S, T ), where S is the curvilinear coordinate along the plate and T is time. A finite number N p of piezoelectric electrodes of length L pi and width H are attached on each side of the plate. For each piezoelectric patch pair, the negative electrodes are shunted across the plate while their positive electrodes are connected to an external circuit. The charge (per unit length in the streamwise direction) in each piezo pair is given by [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF]:

Q i = χ L pi [θ] S + i S - i + CV i , (1)
with V i the voltage between the positive electrodes of the i th piezo pair whose left and right edges are positioned at S - i and S + i respectively. C = C i /L pi is the equivalent capacitance per unit length of the piezo pair and χ is the mechanical/piezoelectrical conversion factor of the piezoelectric material. Considering that each electrode is connected to a purely resistive circuit, Ohm's law leads to:

∂Q i ∂T + GV i = 0, (2)
where G = G i /L pi = 1/R i L pi is the conductivity per unit length of the harvesting circuit, with R i the circuit's resistance.

The voltage V i between the positive electrodes also generates an internal torque in the piezoelectric patches, and thus on the flexible plate, so that the total internal torque in the plate and piezoelectric patches assembly is obtained as

M = B ∂θ ∂S -χ Np i=1 V i F i , (3)
where B is the flexural rigidity, which depends on the respective thickness, Young's modulus and Poisson's ratio of the two constitutive materials [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF]. In equation (3), F i (S) is the polarization function of the i th patch. In the present approach, F i (S) = H s (S -S - i ) -H s (S -S + i ), with H s the Heaviside step function. The conservation of momentum then leads to [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]:

µ ∂ 2 X ∂T 2 = ∂ ∂S   F t τ - ∂ ∂S   B ∂θ ∂S -χ Np i=1 V i F i   n   -P n, (4)
where X(S, T ) is the local plate position and µ is the mass per unit length of the assembly. F t is the internal tension in the plate that enforces inextensibility and P is the pressure force density exerted by the surrounding fluid over the plate. In Eq. ( 4), n and τ are the unit normal and tangent vectors to the plate (see Figure 1). Following previous works on flapping flag modeling [START_REF] Singh | The effect of non-uniform damping on flutter in axial flow and energy harvesting strategies[END_REF][START_REF] Eloy | The origin of hysteresis in the flag instability[END_REF], P is decomposed into a reactive part P reac accounting for potential flow effects (e.g. added mass), corrected by a resistive part P res accounting for the form drag on the plate resulting from lateral flow separation:

P = P reac + P res = m a ρH 2 ẇ -(wu) + 1 2 w 2 κ + 1 2 ρC d H|w|w, (5)
where m a is the added mass coefficient of the transverse section (m a = π/4 for a rectangular plate [START_REF] Blevins | Flow-induced vibration[END_REF]), and C d = 1.8 the drag coefficient [START_REF] Buchak | The clapping book: wind-driven oscillations in a stack of elastic sheets[END_REF]. w and u are respectively the normal and longitudinal components of the plate's velocity relative to the fluid flow, such that Ẋ -U ∞ e x = uτ + wn, with e x the unitary vector in the x direction.

In the present approach, it is assumed that the electrode does not affect the density and bending rigidity of the sandwich plate that are considered homogeneous along the assembly (see [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients[END_REF] for a study on the passive damping of vibrating beams including this effect). Also, in the approach followed in this work, P res accounts for the effect of the lateral flow detachment in the form of a drag. Flow detachment from the trailing edge and wake dynamics are only accounted through the reactive term which models the advection of added fluid momentum (see Ref. [START_REF] Lighthill | Large-amplitude elongated-body theory of fish locomotion[END_REF] for more details). Other models exist that consider the influence of an unsteady wake [START_REF] Alben | Simulating the dynamics of flexible bodies and vortex sheets[END_REF][START_REF] Howell | Interaction between a cantilevered-free flexible plate and ideal flow[END_REF].

A. Non-dimensional equations

Using L, L/U ∞ , U ∞ µ/c and U ∞ √ µc as characteristic length, time, voltage and charge respectively, equations

(1)-( 4) become in non-dimensional form:

q i = v i + α U * γ i [θ] s + i s - i , (6)
β ∂q i ∂t + v i = 0, (7)
∂ 2 x ∂t 2 = ∂ ∂s (f τ τ ) - 1 U * 2 ∂ 2 ∂s 2 ∂θ ∂s n -(M * p res + M * H * m a p reac ) n + α U * i v i δ (s -s - i ) -δ (s -s + i ) n, (8)
with clamped-free boundary conditions. The problem is characterized by the following non-dimensional parameters:

-the non dimensional velocity

U * = LU ∞ µ B , (9)
-the electro-mechanical coupling factor α

= χ √ BC , (10)
-the fluid-solid inertial ratio

M * = ρLH µ , (11)
-the non dimensional plate span

H * = H L , (12)
-the non dimensional electrode length

γ i = L pi L , (13)
-the tuning coefficient of the fluid-solid and electric systems

β = U ∞ C LG . (14)
B. Weakly non-linear form of the equations Equation ( 8) is projected onto the x and y directions in order to obtain two equations for x(s, t) and y(s, t) respectively. The horizontal projection is used to eliminate the tension term f τ from the y projection. Finally, x and its derivatives are eliminated using the inextensibility condition. Keeping terms up to O(y 3 ) one obtains a weakly non-linear equation for y(s, t) (see also Ref. [START_REF] Eloy | The origin of hysteresis in the flag instability[END_REF] for the case of a non-piezoelectric plate):

L(y) + f m (y) + 1 U * 2 f B (y) - α U * f χ (y, v) + M * f res (y) + M * H * m a f reac (y) = 0, (15)
where

L(y) = ÿ + 1 U * 2 y (4) - α U * i v i (δ (s -s - i ) -δ (s -s + i )) + M * H * (y + 2 ẏ + ÿ), (16)
f m (y) = y s 0 ( ẏ 2 + y ÿ )ds -y 1 s s 0 ( ẏ 2 + y ÿ )dsds , (17)
f B (y) = 4y y y + y 2 y (4) + y 3 , (18)
f χ (y, v 1 ..v N p ) = 1 2 y 2 i v i (δ (s -s - i ) -δ (s -s + i )) + y y i v i (δ(s -s - i ) -δ(s -s + i )), (19)
f res (y) = 1 2 C D |y + ẏ|(y + ẏ), (20)
f reac (y) = - 1 2 y y 2 + ẏ y 2 -3y y ẏ -2 ẏ y ẏ - 1 2 y ẏ2 + y s 0 ( ẏ 2 + y ÿ )ds + 2(y + ẏ ) s 0 ẏ y ds -y 1 s y (y + 2 ẏ + ÿ)ds. ( 21
)
Equation ( 16) corresponds to the linearized dynamics while Equations ( 17)-( 21) correspond to nonlinearities related to inertia, stiffness, piezoelectric coupling, resistive and reactive flow effects, respectively.

C. Numerical method

A Galerkin decomposition is used to solve Equation ( 15): the vertical displacement y is expanded as a superposition of clamped-free beam eigenmodes φ p (s),

y(s, t) = ∞ p=1 X p (t)φ p (s). (22)
Next, equation ( 15) is projected on the same set of eigenmodes. Details of the Galerkin projection can be found in A. After truncation to N linear modes, the resulting coupled system of equations is integrated numerically using a semi-implicit step-adaptive fourth order Runge-Kutta method.

D. Harvested energy and harvesting efficiency

Once the dynamics of the piezoelectric flag is obtained using the method above and a limit-cycle oscillation is reached, the total energy harvested by the system is computed as the energy dissipated in all the resistive circuits. In dimensional form, this reads as

P e = Np i=1 G i V 2 i , (23)
with . the temporal average over one oscillation period. The harvesting efficiency η = P e /P f is then computed as the fraction of the fluid kinetic energy flux P f through the vertical cross-section occupied by the flag that is actually harvested and dissipated in the output circuits. Using the non-dimensional variables defined above, this can be written as

η = 1 M * βA * Np i=1 γ i v 2 i , (24)
with A * = A/L and A the flapping amplitude of the flag.

III. FLAPPING INSTABILITY : LINEAR ANALYSIS

As a first validation of our model, we investigate the effect of the piezoelectric arrangement of the flag on its stability for a continuous coverage (s - i+1 = s + i ), and a large number of piezoelectric patches. Beyond a critical value U * c of the flow velocity the system becomes unstable to flutter [START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF]. Figure 2 (left) shows the critical velocity U * c as a function of the mass ratio M * for different values of the coupling coefficient and a large number of small piezoelectric patches (N p = 100). Without piezoelectric coupling, U * c is a decreasing function of M * . Cusps in the curve represent a switch in the mode becoming unstable at the lowest flow velocity (see for example [START_REF] Eloy | Flutter of a rectangular plate[END_REF]). These results match that obtained in Ref. [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF] using a continuous model, thus indicating that for a large amount of piezoelectric patches, the present discrete model indeed converges to the continuous framework.

The electro-mechanical coupling also clearly influences the stability threshold of the flag. In the coupled case (α = 0) and for small values of M * , the damping induced by the electro-mechanical coupling increases the critical velocity U * c . This effect is more important as the coupling factor α is increased. However, for high values of M * , the electro-mechanical coupling has a destabilizing effect, which was interpreted in Ref. [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF] in terms of negative energy waves present on the non-dissipative flag [START_REF] Benjamin | The threefold classification of unstable disturbances in flexible surfaces bounding inviscid flows[END_REF][START_REF] Cairns | The role of negative energy waves in some instabilities of parallel flows[END_REF][START_REF] Doaré | Dissipation effect on local and global stability of fluid-conveying pipes[END_REF].

Several previous studies have focused on the limit of the continuous coverage of the plate by infinitesimal patches (N p 1, [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Bisegna | Optimized electric networks for vibration damping of piezoactuated beams[END_REF]). The discrete approach considered here allows us to investigate the convergence of the results with N p and the applicability of the continuous approach to experimental situations where a finite and often small number of patch pairs is considered. Figure 2 (right) indeed shows a rapid convergence of the critical velocity for N p > 12, beyond which the results of the continuous model are recovered [START_REF] Doaré | Piezoelectric coupling in energy-harvesting fluttering flexible plates: linear stability analysis and conversion efficiency[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF][START_REF] Bisegna | Optimized electric networks for vibration damping of piezoactuated beams[END_REF]. Moreover, it also shows that using fewer patches can significantly modify the stability threshold, positively or negatively, which motivates the present study on a finite number of piezoelectric pairs.

The rapid convergence of the results beyond N p ≈ 10 can be understood as follows: for the values of the mass ratio considered, the dominant unstable mode is of low order resulting in typical wavelengths greater than half a flag length. For a large enough number of patches, the patch length is small enough that the entire mode structure is well captured and the discrete nature of the piezoelectric coverage has no influence.

IV. OPTIMIZATION OF THE ELECTRODES' POSITION

During its self-sustained flapping motion, the flag undergoes complex and non-uniform deformations. Electric charge transfer between the patches' electrodes is intimately linked to the local deformation of the flag. An important challenge for maximizing the amount of energy harvested by such a device therefore lies in the correct positioning of the electrodes. When the electromechanical coupling is large, the piezoelectric patches can significantly modify the flag's dynamics resulting in a non-trivial optimization process.

In order to address this question, we first present a set of preliminary experiments exploring the impact of electrode size and position on the global energy harvesting efficiency of the system. Experimental results are then compared to numerical simulations based on the physical model presented in section II. In the last part of this section, a parametric study is performed based on numerical simulations using more complex geometrical configurations, to identify optimal positioning of the electrodes.

A. Experimental study of the role of electrode size and position

The aim of the experiments described below is to obtain some preliminary evidence of the impact of size and position of the electrodes on the harvested energy. For this purpose we have focused on a simple configuration consisting of a flexible piezo-electric plate covered by a single pair of electrodes of variable length and position, immersed in an axial wind flow of speed U ∞ (figure 3-a). Two different experimental configurations are considered (see figure 3-b for a graphical description of the experiments):

-Experiment 1 : the length of the electrode L p is varied while one of its extremities is fixed at the clamped edge of the plate, -Experiment 2 : the length of the electrode is also varied but with one if its extremities fixed to the free edge of the plate. In these experiments, the plates include a center core consisting of a double sided adhesive tape of length L = 12 cm, width H = 2.5 cm and thickness h s = 100 µm. Both sides of the tape are partially covered by a piezoelectric PVDF film of the length L and width H p = 0.3H, with its negative polarity facing inwards. The thickness of the piezoelectric films is h p = 40 µm for experiment 1 and h p = 50 µm for experiment 2. The remaining free surface at both sides of the scotch tape is covered by a layer of Innova biaxially oriented polypropylene (BOPP) film of exactly the same thickness as that of the piezoelectric film. Figure 3-c shows a schematic view of the plate composition. The piezoelectric PVDF films are covered with an external Cr/Au which serves as the electrode. In both experiments this Cr/Au layer is sectioned in order to adjust the desired value of δ and L p . The objective of covering only one third of the width of the plate with piezoelectric film is to reduce the impact of the sectioning process on the global rigidity of the plate.

For each experiment and a given electrode length and position, the experimental procedure is as follows : the piezoelectric plate is placed in a wind tunnel of rectangular transversal section 10 cm wide and 5 cm high with transparent walls that allow visual access from the outside [START_REF] Doaré | Effect of spanwise confinement on flag flutter: Experimental measurements[END_REF]. The negative electrodes of the piezoelectric films are shunted while a resistance R is connected between the positive electrodes (see figure 3-d). Once the system is set we proceed to generate a wind flow inside the tunnel. Initially at rest, the wind velocity is gradually increased until the plate starts to flap in a self sustained oscillatory regime (upper left image in figure 4), at which point the voltage V o generated along the resistance R is measured (figure 4, upper right). The power dissipated in the resistance is then calculated as the temporal average P e =< V 2 o /R >, for a 5 second time frame. The resistance R is varied in order to find the value of β that maximizes the dissipated power. Indeed, it is expected that the harvesting efficiency or damping presents a maximum between the extreme values β = 0 (short circuit condition) and β = ∞ (open circuit condition). The value of β maximizing the efficiency corresponds to a tuning of the flapping frequency to the characteristic timescale of the RC circuit consisting of the piezoelectric capacitance and output resistor [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. Figure 4 (bottom) shows a good agreement between experiments and simulations obtained from the model described in section II. The discrepancies observed between the experimental and numerical results are mainly due to small differences in the flapping frequency of the flag, mainly caused by the error in the estimated value of the plate rigidity B used in the simulations. It is important to remark that the optimal β value remains fairly constant regardless of the size and position of the electrodes, a direct consequence of the flapping frequency remaining unchanged when changing the electrode position and size. This result was also confirmed in our numerical simulations.

Once the optimal value of β is identified, synchronized voltage measurements and image capturing of the flapping plate are carried out using an oscilloscope and a high speed camera placed over the wind tunnel. Finally, the efficiency of the system is computed as P e /P f (see section II).

Figure 5 shows the evolution of the harvesting efficiency of the piezoelectric flag with the electrode's length L p for both experiments 1 and 2. A good qualitative agreement is observed between the experimental results and the numerical simulations corresponding to the same non-dimensional parameters. For experiment 1 (left panel), the maximum efficiency is achieved when the electrode occupies the whole length of the plate. For experiment 2 (right panel), the maximum efficiency is achieved when the upstream extremity of the patches is located at s ≈ 0.1. Quantitative discrepancies between experimental and simulation results are, as mentioned previously, due to errors in the estimation of the plate rigidity B used in our simulation. Although an effort has been made during the experiments in order to reduce the impact of the electrode sectioning on the plate rigidity, this effect exists, and it is not accounted for in the theoretical model of the plate. Discrepancies in the absolute efficiency values between the two experiments have their origin in the difference between coupling coefficients of each experiment: α = 0.073 for experiment 1 and α = 0.062 for experiment 2.

The good agreement between experiments and simulations also validates the model presented in section II and motivates its use in the parametric study of more complex geometries presented in the remainder of this section.

B. Parametric study

In this subsection, simulation results for the energy harvesting efficiency are presented in the three cases sketched in figure 6:

: I: One single electrode with variable length L p and variable position along the flag surface (two-parameter problem)

: II: Two electrodes of variable length covering the whole surface of the flag (one-parameter problem) : III: Three electrodes of variable length covering the whole surface of the flag (two-parameter problem)

For each case, we analyze the influence of size and distribution of the electrodes on the energy harvesting efficiency, defined as in the experimental case and expressed in non-dimensional form in Equation [START_REF] Lefeuvre | A comparison between several vibration-powered piezoelectric generators for standalone systems[END_REF]. A strong dependence of the efficiency on the system's physical parameters such as the mass ratio M * , β, U * and α is observed [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. Before exhibiting results of geometry optimizations, it is necessary to select a representative set of these parameters. First of all, due to the weakly-nonlinear nature of the model used in the present work, the accuracy of the results is restricted to flow velocities U * near the flapping instability threshold. Far from threshold this model exhibits discrepancies with previous simulations and a fully-nonlinear model of the system should be considered [START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. Thus, all the results presented further in this work correspond to near threshold situations (U * is close to U * c ), and hence, depend on the value of M * selected.

Secondly, two typical values of the mass ratio M * will be considered: M * = 0.6, which corresponds to a plate in a low density fluid (i.e. flag in air), and M * = 10 which corresponds to a plate in a dense fluid (i.e. light plate in water). For each case, the flow velocity U * is adjusted to a typical value just above the instability threshold: U * = 13 and 8.25, respectively.

As explained in the previous section, maximum energy harvesting is achieved when the flapping timescale and the RC circuit timescale compare. Hence, the optimal value of β will depend on the flapping frequency, which in turn depends on the mass ratio M * [START_REF] Guo | Stability of Rectangular Plates With Free Side-Edges in Two-Dimensional Inviscid Channel Flow[END_REF]. On figure 7 the efficiency η in a typical single electrode configuration is plotted as function of β for two different values of M * . One then observes that higher values of M * tend to reduce the optimal value of β. Indeed, while for M = 0.6 the optimal tuning parameter is β ∼ 1 (as in the experimental case), for M * = 10 the value of β is about 0.3. In the following, these optimal values for β are retained.

The influence of the coupling coefficient α, the electrode position and geometry on the flapping frequency has then been tested. No significant impact on the flapping frequency, and consequently on the optimal value of β has been observed. This is expected to result from the relatively low coupling coefficients α considered.

Case I : single electrode

We consider here the case of one electrode of length L p that partially covers the surface of the plate. It is attached at a distance δ from the clamped end. Figure 8 presents the efficiency of the harvesting system as a function of δ and L p , for fixed values of M * = 0.67, U * = 13 and β = 1, and increasing values of the electro-mechanical coupling factor α. Each colored rectangle of these plots corresponds to one single numerical simulation. We observe that regardless of the value of α, the maximum efficiency is found for values of δ = 0 and electrode lengths 0.35 < L p < 0.5. In these cases, the electrode is mainly located on the downstream half of the clamped-free plate. Moreover, the maximum efficiency region is observed to change as the electro-mechanical coupling factor α is increased. For small values of α (∼ 0.1), the maximum efficiency region on the (L p , δ) space extends from long electrodes L p = 0.8 placed at δ = 0.2 to smaller electrodes of length L p ∼ 0.4 and δ ∼ 0.45 (Figure 8a). As the coupling factor α is increased the maximum efficiency is also increased. Simultaneously, the maximum efficiency region shrinks and migrates in the (L p , δ) plane to a new optimal configuration for an electrode of length L p ∼ 0.4 and δ ∼ 0.55 (Figure 8c).

To explain this difference in the optimal configuration, the influence of the geometrical parameters (δ, L p ) on the flapping amplitude is represented on figure 9(a,b) for small and large values of α. These figures show the normalized relative flapping amplitude ∆A N = (A -A 0 )/A 0 in the (δ, L p ) plane for α = 0.1 and α = 0.5, respectively, where A 0 is the flapping amplitude obtained with no coupling (α = 0). Comparison of these figures show that α only impacts the general magnitude of the amplitude fluctuation, and does not modify the parameter values leading to maximum (or minimum) amplitude. In particular, the maximum flapping amplitude corresponds to the optimal efficiency region achieved for high coupling. At large coupling, maximizing the flapping amplitude seems therefore to be the optimal strategy. These maximum flapping amplitudes are also associated with the largest reduction of the critical flow velocity (Figure 9c) which is consistent with the generally observed result that the flag flapping amplitude is a monotonically increasing function of U * -U * c in the vicinity of the threshold. On the other hand, when the coupling coefficient is small, the flapping flag dynamics is only marginally modified by the energy transfers to the piezoelectric patches. In that case, and neglecting the feedback coupling of the piezoelectric, the amplitude of the equivalent current generator (Figure 1a) can be directly obtained using Equation ( 6) from the kinematics of the flag in the case of no coupling. More specifically it is proportional to the relative rotation rate [ θ]

s + i s - i
of the two edges of the piezoelectric patch. Maximum efficiency is therefore achieved for low coupling when ∆θ i is maximum, where

∆θ i = | θ(s + i ) -θ(s - i )| , (25)
a result consistent with Figures 8(a) and 9(d).

Before presenting some results at a higher value of the mass ratio, let us summarize the results obtained so far:

-For weakly coupled systems (α = 0.1), the piezoelectric coupling has no influence on the flapping amplitude and the efficiency scales with ∆θ i (figure 9d).

-For strongly coupled systems (α = 0.5), the piezoelectric coupling α has an influence on the flapping amplitude which in turn also influences the efficiency. It was finally shown that the effect of α on ∆A N can be directly deduced from that of α on ∆U N .

For M * = 10 and α = 0.25, we observe different local maxima of the efficiency in the (δ,L p ) plane. They correspond to an electrode of length L p ∼ 0.2 at δ ∼ 0.3, 0.5 (absolute maximum) and 0.75 respectively. The presence of several local maxima with shorter electrodes is explained by the fact that a larger M * leads to a higher order for the mode selected by the instability [START_REF] Eloy | Aeroelastic instabilities of cantilevered flexible plates in uniform flow[END_REF]: each local maximum corresponds to the localization of the piezoelectric patches in the zones of maximum curvature. This is confirmed by the strong correlation of the maximum efficiency with ∆θ i during the oscillating cycle (Figure 10c).

It should also be pointed out that the maximum efficiency for M * = 10 is about 2.5 times bigger than that for M * = 0.67. This result stems from the difference in mode structures at higher mass ratio: a shorter wave length penalizes long electrodes that only react to the average flag curvature along their length. It also suggests that for higher M * a larger number of smaller electrodes occupying the remaining free space on the plate surface could eventually increase the total efficiency even more. This aspect is investigated further in the following.

Case II: two electrodes

The case of two electrodes covering the entire plate is now considered for the same values of M * , U * and β as in previous case, and depends on a single parameter namely the length of the first piezoelectric electrode x c (see figure 6). For M * = 0.67, the maximum efficiency is obtained for x c ∼ 0.35 (Figure 11 a) and is approximately twice as large as the efficiency obtained for a single electrode covering the entire flag (x c = 0 or x c = 1), and represents a 6% increase of the efficiency in comparison with the optimized one-electrode configuration detailed in Case I. In this new optimal situation, the second electrode covers the most significant part of the plate's surface. Also, it is observed that for the same electrode's length, the downstream electrode leads to the largest energy harvesting, a result that is somehow expected because this region corresponds to the largest deflections of the flag.

For M * = 10, two local maxima are observed around at x c ∼ 0.5 and x c ∼ 0.7. Figure 11(b) confirms that the downstream piezoelectric pair contributes the most to the efficiency for each maximum. However, the maximum efficiency of the two electrodes configuration is slightly smaller than that of case I, a consequence of the full coverage of the plate by the two electrodes, in particular near the leading edge where the local rotation vanishes. Due to the dominant mode structure at large M * , shorter electrodes (either when using three or more electrodes or by choosing to cover only a fraction of the plate's surface) are therefore expected to improve again the efficiency.

Case III: three electrodes

We finally turn to the case of full coverage by three electrodes, described by the positions of the left edge of the first and second piezoelectric patches, x c1 and x c2 respectively. For M * = 0.67 the maximum efficiency is reached for x c1 ∼ 0.4 and x c2 ∼ 0.9 (left plot in figure 12). In this configuration, the second electrode is approximately the same as in the maximum efficiency configuration of case I. The presence of two additional electrodes results however in a 20% increase of the efficiency in comparison with case I.

For M * = 10, the arrangement of maximum efficiency (x c1 = 0.5, x c2 = 0.7) is again related to the optimal configuration of case I, with the second electrode of length L p = 0.2 positioned at δ = 0.5. Taking advantage of the deformation mode, the third electrode of length L p = 0.3 and positioned at δ = 0.7 contributes in almost the same amount to the efficiency as the second electrode. As a result, the total efficiency of the three electrode configuration is about three times higher than that of cases I and II. This clearly demonstrates that for higher mass ratios, a larger number of electrodes results in significant improvements of the efficiency, while this strategy only leads to moderate improvements for lower M * .

V. CONCLUSIONS

This article investigated the possibility of converting kinetic energy of a flow into electrical energy through elastic deformation of a fluttering plate using piezoelectric effect, a problem that has received growing attention in recent years. A critical question lies in the optimization of the energy transfers between the fluid, solid and electric circuit. The present work focused on the role and influence of the location of the piezoelectric electrodes on the efficiency of the harvester, explicitly describing the discrete nature of the piezoelectric coverage of the flag and hence, allowing for a better representation of the experimental and practical situations.

To address this problem, the present work used a combination of theoretical analysis and experiments. A weakly nonlinear model considering a plate with a discrete number of piezoelectric patches pairs was developed and confronted to the experimental results. It was then used to determine geometries that maximize the energy harvesting efficiency in three different test cases: (1) one electrode partially covering the plate, (2) two electrodes fully covering the plate and (3) three electrodes fully covering the plate. For each configuration, the optimal geometries were described and analyzed.

The results presented in this work emphasize the critical role of the piezoelectric elements positioning: because the electric charge transfer responsible for the powering of the output circuit is solely given by the relative change of orientation of the flag between both ends of the patch, the placement of this patch must be carefully designed to take full advantage of the non-uniform flag deformation and its specificities, in particular the dependence of the mode shape with the mass ratio. In the case of a single electrode and small values of the mass ratio, the present study in fact shows that a greater amount of energy can be harvested when positioning the electrodes on the downstream half of the flag. At a higher value of the mass ratio, more than one optimal positioning are found, which is a consequence of more spatial oscillations of the unstable mode shape. Using a larger number of electrodes allows to improve the efficiency of the system particularly for higher mass ratios, where the modal structure corresponds to a shorter wavelength and is therefore better suited to shorter electrodes that are able to capture the flag's deformation best.

The work in this paper extends the results obtained in the context of damping or energy harvesting of vibrating plates in still fluid [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients[END_REF][START_REF] Liao | Optimal placement of piezoelectric material on a cantilever beam for maximum piezoelectric damping and power harvesting efficiency[END_REF] to flow-induced instabilities and more than one piezoelectric patch. The present results allow to conclude that the optimal location for piezoelectric energy harvesting along the plate is on the downstream half, and in that regard are reminiscent of the conclusions of Ref. [START_REF] Singh | The effect of non-uniform damping on flutter in axial flow and energy harvesting strategies[END_REF] obtained for a pure damping model of the energy harvesting process. In both cases, the local deformation of the flag is driving the energy harvesting process, either through changes in the local curvature or in the relative orientation at both ends of the piezoelectric patches, and the response of the deformation of the flag to the presence of energy harvesting plays a critical role in determining the optimal location for large coupling coefficient. Both approaches however differ in the relevant measure of flag deformation (either curvature or relative orientation), which also impacts the detail of the optimal location for the damping or piezoelectric patches.

In the case of weakly coupled systems, it was also shown that the amount of harvested energy can be directly predicted from the kinematics of the flapping flag obtained for zero piezoelectric coupling. Maximizing the harvested energy hence requires to determine the positioning that will ensure a maximum relative variation of the orientation between the two ends of the piezoelectric components. In the case of strongly coupled system, the positioning of the electrodes may influence the flapping amplitude, and consequently the efficiency, because of the feedback effect of the piezoelectric patches on the flag dynamics. Near the stability threshold, the flapping amplitude can be directly related to the critical velocity fluctuations: in that regard, linear stability analysis may be a great tool to estimate the impact of the piezoelectric coverage on the efficiency in the case of large coupling.

FIG. 1 .

 1 FIG. 1. (a) Two dimensional flapping of a flexible plate covered by pairs of piezoelectric patches. (b) Flexible plate flapping in uniform axial flow.

5 M * = 2 M * = 11 M * = 17 FIG. 2 .

 5211172 FIG. 2. (a) : Evolution of the stability threshold as a function of M * for different values of the coupling factor α, β = 1 = cst, Np = 100. (b) : Impact of the number of electrodes on the critical stability threshold. Values of U * c are normalized using the critical value corresponding to Np = 100.

FIG. 3 .

 3 FIG. 3. (a) Description of the experimental setup. (b) Graphic representation of the electrode size and distribution for experiments 1 and 2. (c) Schematic of the multilayer composition of the piezoelectric flag. (d) Connection diagram of the harvesting circuit.

FIG. 4 .

 4 FIG. 4. Top left: Image sequence of the flag deformation over one oscillation period. Top right : Voltage signal generated by the piezoelectric patches during a typical experiment. Bottom : Normalized dissipated power PN as a function of the tuning parameter β. Data corresponds to experiment 1, different symbols correspond to different values of Lp. The normalized power PN is calculated as PN = P/Pmax for each given value of Lp. The solid line corresponds to a numerical simulation using Lp = 1.

2 FIG. 5 .

 25 FIG. 5. (a) : Harvesting efficiency as a function of the electrode length Lp for experiment 1. Experimental parameters are : U * = 14.6, M * = 0.68, H * = 0.3, β = 0.8 and α = 0.073. (b) : Harvesting efficiency as a function of the electrode length Lp for experiment 2. Experimental parameters are : U * = 15.7, M * = 0.61, H * = 0.3, β = 0.8 and α = 0.062. Error bars correspond to the error propagation of the standard deviation of Vo used to calculate the dissipated power Pe =< V 2 o /R > during a 5 second time frame. Solid line corresponds to numerical simulations.

FIG. 6 .FIG. 7 .

 67 FIG.6. Three cases of study in the energy harvesting efficiency of the piezoelectric flag.

FIG. 8 .

 8 FIG. 8. Colormaps of the efficiency η in the (Lp, δ) plane; for all cases, M * = 0.67, U * =13 and β = 1.

9 .

 9 Normalized flapping amplitude difference ∆AN = (A -A0)/A0 in the (Lp, δ) plane for (a) α = 0.1 and (b) α = 0.5 respectively, where A0 corresponds to the flapping amplitude with no coupling (α = 0). (c) Normalized critical fluid velocity ∆U * N = (U * c -U * c0 )/U * c0 as a function of the electrode position in the (Lp, δ) plane where Uc0 corresponds to the flapping critical flow velocity with no coupling (α = 0). (d) ∆θi integrated over one oscillation cycle for a plate with no piezo-electrical coupling. The mode shape of the fluttering flag over one oscillation cycle is depicted in the inset. Other parameters are M * = 0.67, U * =13 and β = 1.

  ∆θ i (α 0) 10. (a) Evolution of the efficiency η with (Lp, δ) for α = 0.25, M * = 10, U * = 8.25 and H * = 0.5. (b) Normalized critical fluid velocity ∆U * N = (U * c -U * c0 )/U * c . (c) ∆θi integrated over one oscillation cycle for a plate with no piezo-electrical coupling. The mode shape of the fluttering flag over one oscillation cycle is depicted in the inset.

* = 10 FIG. 11 .FIG. 12 .

 101112 FIG. 11. Energy harvesting efficiency for a piezoelectric plate covered by two electrodes of variables size for (a) M * = 0.67 and U * = 13, and (b) M * = 10 and U * = 8.25. In both cases α = 0.25 and H * = 0.5. Each figure shows the total efficiency of the two electrode configuration (solid), as well as the efficiency of electrode 1 (dashed) and electrode 2 (dash-dotted).

ACKOWLEDGMENTS

This work was supported by the "Laboratoire d'Excellence" LASIPS (project PIEZOFLAG), and the French National Research Agency ANR (Grant ANR-2012-JS09-0017). S.M. also acknowledges the support of a Marie Curie International Reintegration Grant within the 7th European Framework Program (Grant PIRG08-GA-2010-276762).

Appendix A

The vertical displacement y is expanded on clamped-free beam eigenmodes φ p (s), y = N i=1 φ i (s)q i (t).

(A1)

After truncation to N linear modes, equation ( 15) is projected along mode p. Projection of the forces L(y), f m , f B , f χ , f res and f reac are given in this appendix. Projection of L(y) along mode p reads,

where

For f m (y) we have:

(q i qj qk + q i q j qk )(c p i,j,k + d p i,j,k ), (A5) with

The term f B (y) projects as:

q i q j q k (e p i,j,k + 4f p i,j,k + g p i,j,k ), (A8) where

The mechanical-electrical coupling term f χ gives

with

For the resistive term f res we have :

where

|φ i q i + φ i qi |(φ j q j + φ j qj )ds. (A15)

Notice that unlike the other coefficients t p j,k is a function of time. It is due to the presence of the absolute value. This coefficient must then be recalculated at each time step.

Finally, for f reac we have:

ϕ p i,j,k q i qj q k -N i,j,k=1 p i,j,k q i q j qk (A16) with:


Jean-Paul Caltagirone

Stéphane Vincent
email: stephane.vincent@u-pem.fr

A Kinematic Scalar Projection method (KSP) for incompressible flows with variable density

Keywords: Projection methods, Hodge-Helmholtz Decomposition, Navier-Stokes equation, Incompressible flows

A new Scalar Projection method presented for simulating incompressible flows with variable density is proposed. It reverses conventional projection algorithm by computing first the irrotational component of the velocity and then the pressure. The first phase of the projection is purely kinematics. The predicted velocity field is subjected to a discrete Hodge-Helmholtz decomposition. The second phase of upgrade of pressure from the density uses Stokes' theorem to explicitly compute the pressure. If all or part of the boundary conditions are then fixed on the divergence free physical field, the system required to be solved for the scalar potential of velocity becomes a Poisson equation with constant coefficients fitted with Dirichlet conditions.

Introduction

Solving the equations of incompressible fluid mechanics requires to ensure the coupling of the equation of motion and that of the incompressibility constraint. This coupling can be implicit in a kind of "exact" approach by including the constraint in the linear system as for the method of the augmented Lagrangian [START_REF] Fortin | Méthodes de lagrangien augmenté, Application à la résolution numérique de problèmes aux limites[END_REF]. It proves very efficient and robust however it involves the use of efficient preconditionners associated to iterative solvers and large linear systems due to the coupling of all velocity components [START_REF] Vincent | Augmented lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. application to hydroplaning flows interacting with real tire tread patterns[END_REF]. Another class of methods is to split the resolution of the equation of motion for the application of the incompressibility constraint together with the formulation of an equation for pressure [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF]. One important aspect is immediately apparent to the first authors who have developed numerical algorithms around finite volume methods, that is the spatial location of pressure and velocity unknowns. For collocated variables, instabilities appear and interpolations are needed to mitigate and remove these fluctuations of velocity and pressure [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF]. Harlow and Welch in 1965 have introduced the notion of staggered variables [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF] initially for two-phase flows. This strategy, called the marker and cell method, ensures the coupling, direct or not, of the pressure and velocity fields without disturbance. Since many authors have developed time splitting or prediction-correction methods called projection that involve treating the solving of motion equations and their incompressibility constraint sequentially. Many algorithms allow to obtain convergence orders in time ranging from range O(∆t 1/2 ) to O(∆t 2 ). These orders also depend on the boundary conditions imposed. Some reviews can be found on these techniques including the comprehensive of Guermond et al. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The more recent use of these techniques for the simulation of two-phase flows leads to ill-conditioned linear systems especially for strong density contrasts. Indeed solving a Poisson equation with strongly varying coefficients is very costly in terms of number of iterations of iterative solvers such as conjugate gradient. The direct solvers are performing in two-dimensional space but are unusable in three-dimensional simulations for large numbers of degrees of freedom. Some authors address the problem on an algebraic point of view by specific preconditionning or by the resolution of saddle point [START_REF] Bell | A second order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Benzi | Preconditioning techniques for large linear systems : a survey[END_REF][START_REF] Benzi | Numerical solution of saddle point problems[END_REF][START_REF] Cai | Efficient variable-coefficient finite-volume Stokes solvers[END_REF]. Other ways are sought for example by Guermond et al. [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure poisson equation[END_REF] [16] extracting the density of the Poisson equation. This approach is effective mainly for small density ratios. It can be used with some caution for flows involving open boundary conditions [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. The recent fast pressure-correction method of Dodd and Ferrante [START_REF] Dodd | A fast pressure-correction method for incompressible two-fluid flows[END_REF] is also based on the factorization of the density in the projection step by using of the minimum density between two separated fluids with the introduction of a pressure source term in the Poisson equation. This method has been validated against standard capillary test cases and it was utilized to simulate the interaction between a homogeneous isotropic turbulence and 6260 spherical particles. Works based on a vector approach of the resolution of the projection step [START_REF] Caltagirone | A vector projection method for solving the Navier-Stokes equations[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] are particularly effective for the simulation of two-phase flows with large density contrasts. However it requires the solving of large linear systems induced by the coupling of all velocity components in the projection. In the same field, parallel work centered on the discrete Helmholtz-Hodge decomposition bring potential solutions to using it for solving partial differential equations such as Navier-Stokes equation [START_REF] Ph | Fast discrete Helmholtz-Hodge decompositions in bounded domains[END_REF][START_REF] Lemoine | Discrete Helmholtz-Hodge decomposition on polyhedral meshes using compatible discrete operators[END_REF]. Other potential applications of this decomposition are numerous. They are detailed in the review of Bhatia et al. [START_REF] Bhatia | The Helmholtz-Hodge decomposition -a survey[END_REF]. The approach proposed here is based entirely on the mechanics of discrete media [START_REF] Caltagirone | Discrete Mechanics[END_REF]. This formulation of the momentum conservation equation results in a set of equations that is very different from the standard Navier-Stokes equations. Its constitution resumes from the fundamental law of dynamics, Newton's second law, and on a vision of differential geometry. Thus a discrete equation of motion is obtained in the form of a natural decomposition of Helmholtz-Hodge in irrotational and solenoidal parts. This approach has led to a purely vectorial version of the projection method where at the end of this second step the boundary conditions of the problem are satisfied [START_REF] Caltagirone | Application de la décomposition de Hodge-Helmholtz discrète aux écoulements incompressibles[END_REF]. This short note is devoted to a scalar version of vectorial projection where irrotational component of the velocity is sought by a Helmholtz-Hodge decomposition of the scalar potential of velocity. A detailed presentation of the projection algorithm is first given and several illustrative examples of flows with varying densities are provided for discussion.

Kinematic scalar projection (KSP) method

The resolution of the equation of motion in an incompressible formulation associated with the boundary conditions of the physical problem is the objective of projection methodologies. This motion equations can be the Navier-Stokes equations or the equations coming from the discrete mechanics [START_REF] Caltagirone | Discrete Mechanics[END_REF]. For flows at constant density and at constant viscosity both formulations are equivalent whereas it is not the case for variable fluid properties encountered in two-phase flows for example. As this choice is regardless for the scalar projection method under consideration here, the Navier-Stokes equations of motion for an incompressible flow are resumed

     ρ dV dt = -∇p -∇ • µ ∇V + ∇ t V ∇ • V = 0 (1)
where V is the velocity, t the time, p the pressure, ρ the density and µ the dynamic viscosity.

The boundary of the physical domain Ω is noted Σ. By decomposing the particular derivative of velocity and discretizing equations at first order in time (a second order is easily obtained with the present method by using a second order Taylor expansion in time), a prediction step can be formulated for the intermediate velocity

V *      ρ V * -V n ∆t -ρ V n • ∇V * = -∇p n -∇ • µ ∇V * + ∇ t V * V * × n = V n+1 ΨΣ × n on Σ (2)
with p n the pressure at time n∆t, ∆t the time discretization index and V Ψ the solenoidal component of the desired velocity, that is V n+1 for a time discretization for which t n+1 represents the time at the end of the prediction and projection steps of the time splitting approach. Twophase flows with variable density encompass a wide variety of different physical situations. For example, the liquid-particle gas flows do not present the same difficulties as the hydraulic jump flows even if the density ratio is the same. In the second case, this is the difference in density associated with the gravity that generates the flow motions. In many cases, the algorithms described in the introduction section are sufficient to simulate the physical behavior of the problem. The Helmholtz-Hodge decomposition of the predicted velocity field is sufficient to maintain the balance between the effects of gravity and the dynamic effects even if the field V Φ is not exactly an irrotational field. This is the case for example for the natural convection presented in section 3.2. In the absence of the incompressibility constraint in the equation of motion in the prediction step, the divergence of V * is not zero and only the normal velocity component is respected at the boundaries. The predicted field V * = V Φ + V Ψ includes both solenoidal and irrotational contributions. The divergence free component V ΨΣ is known on Σ thanks to the physical boundary condition to impose. However, the projection step does not allow to maintain it as it is related to a scalar equation. The scalar projection step consists in searching the irrotational component V Φ = ∇Φ and then obtain the divergence free component by difference with V * . One of the Helmholtz-Hodge decomposition methods of the velocity V * = ∇Φ + ∇ × Ψ into its two divergence and rotational free components amounts to applying the divergence operator to V * and ∇Φ that are equal since the divergence of the rotational is zero. It can be demonstrated [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] that null flux conditions on the scalar potential have to be associated to the Poisson equation :

   ∇ 2 Φ = ∇ • V * ∇Φ • n = 0 on Σ (3)
for Dirichlet boundary conditions on V * . Other types of boundary conditions can be also considered as discussed for example in [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF] [14] [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. In fact it can be built as many irrotational fields as boundary conditions applied to the system. However, there is little chance of finding a solenoidal field that satisfies the boundary conditions of the physical problem. This step is purely kinematics and does not result from any numerical time splitting. It does not involve either the density which is perfectly legitimate. The flow may be at continuously variable density or two-phase, there is no physical reason or mathematical argument that leads to associate the density in the scalar potential velocity Φ. The dynamic part of the projection step contrariwise brings up the density for the pressure p. The following of the section specifically returns to the link between the scalar potential of velocity and pressure.

The first equation of system ( 2) is a prediction step. It can be solved with the physical boundary conditions of the problem. Its solution V * does not satisfies the divergence free constraint. Indeed, the solving of this prediction step with a scalar potential p n not adapted to the boundary conditions of incompressibility and unknown at the solving time of the prediction step introduces a non zero irrotational component that has to be removed by means of a Hodge-Helmholtz decomposition. The second equation of ( 1) cannot be solved directly as the actual component V Φ of the actual field has to be deduced from the difference between the fields V * and V n+1 Ψ obtained after the decomposition of the field V * . It can be observed that

V Φ = V n+1 Φ derives from a potential V Φ = ∇Φ.
The second equation of ( 1) can in this way be directly and explicitly integrated [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]. The scalar potential of the velocity is obtained by a projection on a divergence free field by solving (3). If we discard the inertial and viscous terms and consider the momentum equation expected to be solved and the prediction step really solved, we have

         ρ V n+1 -V n ∆t = -∇p n+1 ρ V * -V n ∆t = -∇p n (4)
By substraction of these two equations, we obtain

ρ V n+1 -V * ∆t = -∇ p n+1 -p n (5)
We also have

   V * = ∇Φ + ∇ × Ψ = V Φ + V Ψ = V Φ + V n+1 ∇ • V n+1 = ∇ • (V * -V Φ ) = ∇ • (V * -∇Φ) = ∇ • V Ψ = 0 (6)
By combining ( 5) and ( 6), an equation linking the pressure increment p ′ = p n+1 -p n and the irrotational velocity component can be written

ρ V Φ ∆t = ∇p ′ (7)
By using ( 3) and ( 7), we finally get the KSP projection step

             ∇ 2 Φ = ∇ • V * V Φ = ∇Φ ρ V Φ ∆t = ∇p ′ (8)
The pressure increment is obtained by considering the Stokes theorem

b a ∇Φ • tdl = Φ b -Φ a (9)
According to ( 8) and ( 9),

p ′ b -p ′ a = b a ∇p ′ • t dl = b a ρ V Φ ∆t • t dl (10)
Equation ( 10) is correct only if ρ V Φ ∆t is a gradient. On a discrete point of view, ρ is constant over each mesh edge and ∆t is also constant, so that ρ

V Φ ∆t = ρ ∇Φ ∆t = ∇ ρΦ ∆t is clearly a gradient
and expression [START_REF] Caltagirone | A vector projection method for solving the Navier-Stokes equations[END_REF] always hold even for multi-phase flows. Finally, the update of the velocity and pressure fields are given by expressions :

       V n+1 = V n+1 Ψ = V * -V Φ p ′ n+1 b = p ′ n+1 a + b a ρ ∆t V Φ • t dl (11)
where a and b correspond to the edges of the segments Γ forming the computational mesh where the density is constant. This last relation corresponds to the application of one of the forms of the Stokes theorem which allows to update the potential geometrically if the velocity field is irrotational what is. In the case where the segment Γ is intersected at point c by an interface separating fluids of densities ρ 1 and ρ 2 , the Stokes formula becomes :

p ′ n+1 b = p ′ n+1 a + c a ρ 1 ∆t V Φ • t dl + b c ρ 2 ∆t V Φ • t dl (12)
As for the velocity, it stays continuous and constant along all the segment. Point c will be determined thanks to an interface tracking method of VOF, Front-Tracking or Level-Set type. The pressure at time n + 1 is then p n+1 = p n + p ′ . Practically, we start integrating the pressure increment from an arbitrary mesh point for which ρ is known (belonging to one phase or another) and by stating p ′ = 0 at this point. This procedure is valid on structured or unstructured grids. The solving algorithm finally obtained is very simple to implement :

-1 -Prediction step : solving of the first equation of system (2) to obtain V * -2 -Projection step : decomposition of the field V * in a gradient of the potential Φ and a rotational part by using system (3) -3 -Estimate of the irrotational component V Φ = ∇Φ and of the solenoidal component by using the difference

V n+1 Ψ = V * -V Φ
-4 -Update of the velocity and pressure by considering [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF].

The solution from one iteration at time t n+1 consists of the pressure and divergence free velocity. It does not meet the imposed tangential physical boundary conditions as do the other projection methods. A boundary layer is so created that disappears during the time iterations with the imposition of the physical boundary conditions in the two stages of prediction and correction. The thermodynamic pressure p is in fact utilized only to evaluate the properties of the fluids such as the density and the viscosity. For multi-phase flows at constant densities, the Bernoulli pressure is well adapted in this case. The present KSP algorithm, also called DSP by [START_REF] Caltagirone | Application de la décomposition de Hodge-Helmholtz discrète aux écoulements incompressibles[END_REF] in the framework of discrete mechanics equations, inverts the calculation steps for and velocity compared to standard projection approaches. For classical projections methods, the pressure is first estimated as the solution of a Poisson equation with variable coefficients and the velocity is then explicitly obtained. In the KSP method, the velocity is first decomposed by a purely kinematic processus and the pressure is then updated by the explicit application of the Stokes theorem. The solving of a Poisson equation with variable coefficients is a difficult task whose complexity increases with density ratios, especially with large grids on massively parallel computers, whereas the KSP method is not sensible to these density variations.

When high density gradients are associated to large magnitude source terms, it can be necessary to perform a preliminary Helmhotz-Hodge decomposition of source terms s, before the time evolution loop begins. In this way, the initial condition for pressure at mechanical equilibrium is then built as

-∇p o + s Φ = -∇p o + ∇Φ s = 0 (13)
The initial condition is given by

p o = Φ s . The source term s = s Φ + s Ψ = ∇Φ s + ∇ × Ψ s .
s Ψ = ∇ × Ψ s does not directly induce a pressure. However, it generates the motion through the increase of the solenoidal velocity field. As a consequence, s Ψ has to be kept in the momentum equations instead of s when initial condition (13) is implemented. On the contrary, a source term deduced from a gradient field does not involve any flow motion. In this configuration, it is directly integrated inside the pressure gradient to form a new pressure field. The decomposition of the initial source term s is achieved by using a Poisson equation similar to (3) so as to obtain Φ s :

   ∇ 2 Φ s = ∇ • s ∇Φ s • n = 0 on Σ (14)
Prior decomposition of the source term eliminates the adverse effects induced by exchanges between the pressure effects and all the other effects (viscous, inertial) that cause unwanted local and instant acceleration that affects the quality of the two-phase behavior. This KSP version suitable for very constrained two-phase flows requires an additional projection step and the computing of the solution of a Poisson equation with constant coefficients. However, it allows to build a very robust algorithm. The efficiency of solvers with constant coefficient Poisson equations offsets the additional cost of the resolution of an additional equation with respect to the conventional methods of projection. In general the introduction of a significant source term in the Navier-Stokes equations generates difficulties due to the destabilization of the vector field by the scalar potential that it contains whereas the latter do not participate to the movement itself.

3 Illustration test cases

Static equilibrium between two fluids under gravity effects

The considered problem is very simple, it consists of a square cavity of unit height filled with two immiscible fluid whose densities are ρ 1 and ρ 2 . It is assumed that the two fluids are initially separated and the heavy fluid (1) occupies the lower half of the cavity. The stationary solution is simple : the velocity V is zero and the pressure field satisfies p(y) = p 0 + ρ(y) g • e y y. The initial pressure field is zero. The walls of the cavity are assumed impermeable and adherent, V Ψ = 0. Details of the different steps of the time splitting algorithm on this problem are the following. In the absence of initial velocity, this velocity field derived from the prediction step (1) is V • e y = -∆t g. If this field is divergence free within the cavity that is not the case near horizontal walls due to boundary conditions. Assuming incompressibility, these variations of divergence restore a linear distribution of the scalar potential of the velocity along the vertical axis. The numerical solution of Φ obtained up to a constant by solving the Poisson equation ( 3) is represented by its evolution along y as

Φ = -∆t g y (15)
The irrotational velocity field is V Φ • e y = ∇Φ • e y = -∆t g. The difference V * -V Φ = 0 is the requested velocity V n+1 = V Ψ = 0. The pressure update by the Stokes theorem is then

p n+1 = p 0 -ρ(y) g y (16)
where p 0 is the reference pressure chosen in an arbitrary manner. At the end of the two stages the theoretical solution is obtained exactly up to computer accuracy. This is the sought two-phase hydrostatic equilibrium. Figure [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] represents the opposite of the evolution of the pressure along y for two density ratios ρ 1 /ρ 2 = 4 and ρ 1 /ρ 2 = 10 6 . In the present case the two fluids have constant densities and each pressure point is in a fluid or in the other, the density is absent from differential operators and appears only for the increase of the pressure. The first phase for determining the potential Φ is independent of the density variations and the projection phase being explicit and local, the solution will always be accurate. All two phase incompressible flows can be simulated with the KSP method with the same efficiency.

Natural convection in a differentially heated cavity

Flows with variable density can be very different in nature, flows involving several immiscible phases, flows with phase changes, etc. Flows with continuously varying density which can be approached in the context of the incompressible approximation belong to this class. Natural convection is an example especially when the temperature differences are important and when the Boussinesq approximation is no longer valid. The example below aims to show that the proposed methodology allows to find precisely the solution adopted by many authors after multiple comparisons. This is the case of a cavity filled with air subjected to a horizontal temperature gradient in a gravity field. Natural convection induced by density variations is quantized by the Rayleigh number and the Prandtl number. The selected configuration correspond to a value of the Rayleigh number such that Ra = 10 5 and it admits a stationary solution. Nusselt number that characterizes the heat transfer between the two isothermal walls is the main result of the problem. The reference solution is obtained by a finite volume method on a Cartesian staggered mesh with augmented Lagrangian technique [START_REF] Vincent | An adaptative augmented lagrangian method for three-dimensional multi-material flows[END_REF] to ensure incompressibility constraint. The spatial order of convergence of the Nusselt being strictly equal to 2, it is possible, using Richardson extrapolation [START_REF] Richardson | The deferred approach to the limit[END_REF], to derive the reference value of Nusselt for N → ∞. The present test case Reference N → ∞ KSP N = 1024 2 Nusselt 4.521638 4.521614

Table 1. Natural convection in a differentially heated cavity for a Raleigh number of 10 5 . The reference Nusselt number is obtained by Richardson extrapolation [START_REF] Richardson | The deferred approach to the limit[END_REF]. The value of the KSP method correspond to solution on a 1024 2 mesh.

is almost trivial but it has the advantage of providing a reference for flows with low variable density in an incompressible formulation. Furthermore, the Nusselt number is very sensitive to the numerical methodology. It allows anyway to find precisely a well-known solution with an method.

Sloshing in a 2D tank

With the addition of specific source terms, the system (1) can model many phenomena according to external actions such as gravity, capillary forces or rotation. In the case of a constant and uniform force of gravity, surface gravity waves of different nature can grow and maintain over large time constants at a fluid/fluid interface. This is the case of solitary waves or swell. In the present test case, a liquid sloshing in a cavity partially filled of gas is considered. First order involved mechanisms are inertia and gravity. Both although formally compressible fluids give rise to a motion that can be considered as incompressible at large time, so that the KSP method can be Consider a cavity of length L and height H that contains a fluid of density ρ 2 and viscosity µ 2 topped with a fluid of density ρ 1 viscosity µ 1 . The interface between the two immiscible phases is slightly disturbed in a sinusoidal manner such that its initial height h(x) is given by h

(x) = H/2 + A cos(πx/L) (17)
with H = 0.1, L = 0.1 and A = H/100 in linear regime and A = H/3 in non-linear regime. Under the effect of gravity, the interface oscillates around an equilibrium position, i.e. an horizontal reference line. At equilibrium, the lower fluid occupies a height H/2. Figure 2 shows the time history of vertical interface position during time in linear regime. The amplitude of the initial perturbation permits to stay within the framework of [START_REF] Landau | Mécanique des Fluides[END_REF]. The notion of viscosity contributes nothing to the problem, only a more or less rapid decay of the waves would be obtained according to the value of viscosities. These are considered null, mu 1 = mu 2 = 0 and the diffusion term of the momentum disappears from the equation of motion. The evolution in time is thus conditioned by competition between the inertia of the fluid determined by the term ρ V • ∇V and gravity source term s = ρ g. The latter term is not derived from a scalar potential and Hodge-Helmholtz decomposition of s = ρ g = ∇Φ s + ∇ × Ψ s highlights two non-zero contributions thereof. The irrotational part modifies the scalar potential p o which includes the static gravitational effects and the vector potential changes the mechanical equilibrium. Coupling with inertia causes the sloshing movement whose frequency may be calculated by the linear theory. If the initial disturbance of the interface is defined by Fourier modes, i the longitudinal direction and j for transverse modes, the linear theory allows to express the frequency [START_REF] Lamb | Hydrodynamics. 6 e edition[END_REF] :

f = 1 2π π g i 2 L 2 + j 2 l 2 th πH i 2 L 2 + j 2 l 2 (18)
where l is the width of the domain along y. In two-dimensions, j = 0 and l = 1. We also define the pulsation ω and period T :

ω = 2 πf = 2π T (19)
The expression of the theoretical frequency [START_REF] Lamb | Hydrodynamics. 6 e edition[END_REF] was established from a linear stability theory for a fluid density ρ 2 in the absence of fluid located above. When the densities ρ 1 and ρ 2 are close, it is necessary to introduce a correction [START_REF] Landau | Fluid Mechanics[END_REF] which gives the relationship Selected fluids are water and air and the corresponding densities are ρ 2 = 1000 kg.m -3 and ρ 1 = 1.1728 kg.m -3 . Only the first 2D mode is tested, i.e. i = 1 and j = 0.

f = 1 2π π g i 2 L 2 + j 2 l 2 ρ 2 -ρ 1 ρ 2 + ρ 1 th πH i 2 L 2 + j 2 l 2 (20)
The time step is equal to 10 -3 s which achieves sufficient accuracy on the frequency of oscillations. Figure 2 shows the periodic changes in the height of the fluid 2 on one edge of the Theory Simulation Period 0.3742 0.3748 Table 2. Sloshing periods in a square cavity for the first linear mode. field x = 0. Gravitational forces introduce a downward movement of the area where the free surface is highest. In the absence of viscous forces, the oscillatory motion is governed by the confrontation between gravity and inertia. It is observed that the oscillations persist for a long time without significant attenuation. It is also to be noted that the wave attenuation is even lower than the time step decreases. This check is used to test the entire methodology : the equation of motion, KDP time splitting algorithm, time and space discretization, interface tracking, etc. To quantify the errors introduced by the different modeling and discretiezation steps, frequency numerically obtained is compared with the theoretical frequency formulated by relation 18. Table 2 rather presents the period of the oscillations. There is a very good agreement between simulations and theory. This test validates the KSP approach in the presence of source terms. This example also serves to show that the formulation conserves kinetic energy when the viscous effects are neglected. Although in this case no transfer of momentum by viscosity is possible, that does not mean that the curl of the velocity field is zero.

The non-linear mode is illustrated in figure 3. It is observed that non-linear interface deformations are nicely handled by the KSP method with vertical interface oscillations being submitted to irregular sinusoidal modes.

Rotation flow

The present test case corresponds to a block in rotational flow in a cylindrical cavity of radius R. The steady rotational velocity Ω o is constant and the tangential component of absolute velocity is V a = v θ e θ = Ω o r e θ . The cavity is filled with two fluids of density ρ 1 and ρ 2 and the interface is initially located at r = R/2. The viscosity has no influence at least in the absence of differential movement relative to the plug flow. The flow motion can be treated in the moving frame relative to Oz axis. In the present configuration, the relative velocity V is chosen equal to zero in the whole domain. According to the momentum equations, the equation for pressure is given by Figure 4. Comparison between the scalar projection (PS) and the kinematic scalar projection (KSP) on a triangular unstructured mesh. From left to right are represented the rotationnal of velocity for classical scalar projection (PS) with ∇ × V ± 0.1, the rotationnal of velocity KSP method with ∇ × V ± 10 -15 , the scalar potential Φ and the pressure p given by the KSP method. The divergence of the velocity is zero to almost computer error for both methods. The inner radius of the circle (black line) is R = 0.5 m and the densities are ρ 1 = 1 kg.m -3 inside the circle and ρ 2 = 4 kg.m -3 outside of it.

ρ Ω × Ω × r = -∇p (21)
with Ω = ∇ × V. The pressure can then be calculated analytically and can then be compared with the numerical solution

p(r) = ρ Ω 2 r 2 2 + p o ( 22
)
where p o is a selected constant chosen equal to zero on the axis. Since the density is not constant in the whole area, the pressure field will be calculated in the two fluid sub-domains on an analytical point of view. The KSP method is now applied from a zero velocity field V = 0 and a zero pressure field p = 0. Equation ( 2) applied to the problem gives the prediction rate V * which is not divergence free. This is a centrifugal velocity oriented outward as

V * = ∆t Ω 2 o r (23)
From that predicted velocity field, it is possible to apply the projection phase (3) for obtaining the scalar potential Φ of the velocity. However, in the present test case V * is a gradient field from the potential Φ defined in a constant, V * = ∇Φ that satisfies

Φ = ∆t Ω 2 o r 2 2 (24)
As the theoretical solution is a polynomial of order two, it is expected that the numerical solution will be accurate. Indeed, all polynomial of order lower or equal to two can be represented exactly by a spatial discretization scheme of order equal to two. Solving the Poisson equation [START_REF] Bell | A second order projection method for the incompressible Navier-Stokes equations[END_REF] actually gives the expected result, as reported in figure 4. As the projection is purely kinematics, the correction velocity V φ = -∇Φ is completely continuous and has no discontinuities at the interface between the two fluids. The numerical solution obtained with KSP at the end of the two stages is zero. The pressure is obtained from the Stokes formula [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] taking care to calculate the integral by piece if the segment Γ of a given mesh cell is intersected by the interface. Here the interface is known analytically and it is simple to specify the position of the intersection point on the segment. Not only the divergence of V n+1 but also the velocity field itself are zero with the KSP method, in agreement with the expected analytical solution. This is not the case with standard projection methods that generate rotation velocity components that are non physical. The solution to this problem was obtained with the KSP algorithm accurately regardless of the mesh type, i.e. structured or unstructured, in one time iteration consisting of a prediction step and a kinematics correcting step.

With classical scalar projection methods, the velocity is calculated from the pressure correction V n+1 = V * -ρ/∆t∇p ′ . For flows with variable density, the divergence of this velocity is zero but its curl is not zero anymore. Indeed, if ∇p ′ is a gradient, ρ /∆t ∇p ′ is not. This is not consistent with what is obtained by using the prediction V * in KSP. In the case where the velocities are large, this residual curl is merged with that of the flow, but for low flow rates or for equilibrium situations like the case shown in this section, it is necessary that V n+1 is also a gradient field (rotation or a drop in capillary effect) in order to satisfy equilibrium state. In all other two-phase cases, the PS numerical technology generates an artefact such as some spurious rotational.

Conclusions and discussions

The kinematic KSP projection method for solving the equation of discrete incompressible fluid motion [START_REF] Caltagirone | Discrete Mechanics[END_REF] essentially solves various problems of incompressible flows, including flows with significant density variations. Unlike conventional methods where the pressure is first calculated from a Poisson equation with variable coefficients, the irrotational velocity is calculated first in KSP. The scalar potential of velocity is then obtained by solving a Poisson problem with constant coefficients that is insensitive to density variations. The scalar potential of the amount of acceleration, i.e. the pressure, is obtained by the Stokes' theorem by introducing at this stage the local density. In terms of accuracy in time and space, the results are very close to those of the conventional projection methodology for flows at variable density. However, the large variations in density introduce local consistency defects in standard projection methods due to interpolation of density at the location of each component of the velocity. The pressure undergoes non physical variations that can lead to unstable or not physical behaviors. The KSP method allows to find consistency between the pressure and the local density. This method can be interpreted as a simple splitting of the motion equation of the continuum mechanics previously discretized in time. It is based on an original formulation of the law for fluid dynamics written as a discrete Helmholtz-Hodge decomposition. The proposed KSP time splitting approach satisfies the following properties :

-the continuous media properties of differential operators, i.e. ∇×∇Φ = 0 and ∇•∇×Ψ = 0, are satisfied to almost computer error -the space convergence order is 2 with a centered scheme and the time convergence order can be 1 or 2 depending on the order of the Taylor expansion used for the time derivative of the momentum conservation equations -the numerical solution is exact whatever the mesh for all theoretical solution of order equal of less than 2 -this is a prediction-correction method whose artefact are well known, i.e. artificial boundary layers are generated by the projection step near the boundaries. Their magnitude decrease during time iterations -Unlike conventional projection methods, the resolution steps for pressure and velocity are reversed. The scalar potential of the velocity Φ is first obtained and then the physical potential, i.e. the pressure, is updated explicitly and accurately -the Poisson equation for velocity potential is at constant coefficients and the velocity potential does not depend on density -the solving of the linear system is easy and allows the use of existing efficient parallel solvers As a conclusion, the KSP method is, among those existing in the literature, the easiest method to implement since it consists in solving a Poisson equation with constant coefficients.

Figure 1 .

 1 Figure 1. Test case of static equilibrium of two fluids of different densities. The opposite to the pressure -p is plotted. In red, ρ 1 /ρ 2 = 4 and in blue ρ 1 /ρ 2 = 10 6 . The spatial approximation is N = 8 grid cells. The solution obtained in one time iteration is exact to almost computer accuracy, i.e. up to 10 -15 relative error.

Figure 2 .

 2 Figure 2. Sloshing of a sinusoidal wave in a 2D tank in linear mode -Time history of vertical interface position at x=0 m.

Figure 3 .

 3 Figure 3. Sloshing of a sinusoidal wave in a 2D tank in non linear mode -initial interface shape (top left), interface solution after s (top right) and time history of vertical interface position for x= 0 m (bottom).