T M Faure

) K Roncin

B Viaud

T Simonet

L Daridon

Flapping wing propulsion: comparison between discrete vortex method and other models

Keywords: bio-inspired concept, multi-physics, aeroelasticity, discrete vortex method, detached flow

The cetaceans propulsion by a periodic flapping motion of their fluke is considered and studied on a benchmark flexible straight wing. The aim of the study is to validate low order models for this configuration. First, the two-dimensional rigid case is investigated, comparing the aerodynamic performance of the airfoil periodic motion versus the reduced frequency, with published data and unsteady Reynolds-averaged numerical simulation results. It appears that viscous drag modeling must be added to the discrete vortex method, in order to obtain sensible thrust results, for Garrick frequencies below 2. All high and low order models agree at the remarkable Garrick frequency of 1.82, although the experiment shows a lower efficiency of about 25%. The positions of the shed vortices match comparing the unsteady Reynolds-averaged numerical simulation and the discrete vortex method. Then, the three-dimensional leading-edgesuction-parameter modulated discrete vortex method is extended, by means of a lifting line theory. A modification of the method is proposed in order to consider wing dihedral, resulting from the spanwise flexibility. The configuration considers a reduced frequency of 1.82. Three types of spanwise wing flexibility are examined. For the inflexible and flexible cases, a reasonable agreement is observed between the different methods for each coefficient. The intermediate flexible wing provides a better thrust coefficient, while excessive flexibility proves to be detrimental. Vorticity fields are compared with previously published data for the three wings. For the highly flexible wing and the right choice of deformation parameters, the discrete vortex method produces reliable results.

I. INTRODUCTION

The periodic vertical plunging of a fluke is a means of propulsion for the cetaceans with efficiencies between 83% and 92% [START_REF] Liu | Propulsive performance from oscillating propulsors with spanwise flexibility[END_REF][START_REF] Fish | Comparative kinematics and hydrodynamics of odontocete cetaceans: morphological and ecological correlates with swimming performance[END_REF][START_REF] Fish | Hydrodynamic performance of aquatic flapping: efficiency of underwater flight in the manta[END_REF][START_REF] Ayancik | Three-dimensional scaling laws of cetacean propulsion characterize the hydrodynamic interplay of flukes' shape and kinematics[END_REF] . For a global overview, testifying of the continuous interest in the subject, see the recent review work by Wu et al. [START_REF] Wu | A review on fluid dynamics of flapping foils[END_REF] . This type of thrust generation is a source of inspiration for the development of innovative propulsion systems, for water or aerial drones, which requires a better understanding of the fluid dynamics involved in the propulsion mechanism. In order to study the flapping performance, for various parameters, such as the wing aspect ratio and its frequency, low computing cost simulation methods are necessary. They need a preliminary validation from published data. That is the reason why a benchmark configuration is studied hereafter.

The oscillatory motion in pure heave of an airfoil, for propulsion purpose, is the subject of many publications. Most of them consider a two-dimensional (2D) flow, as is done in the very first analytical approach for an incompressible and inviscid flow [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] . In that study, the time-averaged thrust coefficient and lift peak coefficient are obtained using a linearized potential flow theory. Lewin & Haj-Hariri 7 investigated the flow features, developing with the oscillation frequency and heave amplitude, in an incompressible viscous flow by a numerical model, where the governing equations are discretized using a conformal map. Downstream of the airfoil, different vortex shedding regimes are identified from the vorticity fields, a) Electronic mail: thierry.faure@ecole-air.fr. depending on the heave amplitude and frequency. Kamrani et al. [START_REF] Kamrani Fard | A leading-edge vortex initiation criteria for large amplitude foil oscillations using a discrete vortex model[END_REF] developed a discrete vortex method, with a leading edge vortex initiation criterion, for the sinusoidal heaving and pitching motion of an airfoil. It is based on a panel method, and an empirical trailing-edge separation correction is implemented. The range of Reynolds numbers considered in the study is 1000 -45, 000 and the reduced frequency is between 0.06 and 0.16. The sinusoidal oscillation of a NACA 0012 airfoil is studied by Young & Lai [START_REF] Young | Oscillation frequency and amplitude effects on the wake of a plunging airfoil[END_REF][START_REF] Young | Numerical simulation of the unsteady aerodynamics of flapping airfoils[END_REF] , for a Reynolds number of 2 × 10 4 , using a 2D Navier-Stokes solver. The development of the thrust and power coefficients and of the propulsive efficiency is given versus the plunging frequency and amplitude. Spanwise flexibility of the wing is considered by Liu [START_REF] Liu | Propulsive performance from oscillating propulsors with spanwise flexibility[END_REF] with an unsteady panel method. An experimental characterization of the flow is made in a water tunnel [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] , for a wing built on a NACA 0012 airfoil with an aspect ratio of 6 and for Reynolds numbers in the range 10, 000 -30, 000. Three wings with different stiffnesses are tested, and a 50% thrust benefit is observed for the wing of intermediate flexibility, while excessive flexibility proves to be detrimental. The very same wings are reported in high-fidelity aeroelastic computations [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] . These authors find the same conclusions regarding the advantage of a moderate flexibility, and produce data of the spanwise development of the flow on the wing upper surface. Recently, a non-linear three-dimensional (3D) lifting-line method is implemented for wings of similar stiffness [START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] . The interaction between the leading edge suction and the shear layer developing on the airfoil during dynamic stall is modeled [START_REF] Deparday | Modeling the interplay between the shear layer and leading edge suction during dynamic stall[END_REF] . Wind tunnel experiments are performed on the effects of passive, inertiainduced surface deformation, at the leading and trailing edges, of an oscillating airfoil energy harvester [START_REF] Siala | Experimental study of inertia-based passive flexibility of a heaving and pitching airfoil operating in the energy harvesting regime[END_REF] .

The present investigation applies the Leading-edgesuction-parameter-modulated Discrete Vortex Method (LDVM) and other models to the oscillatory heave of a propulsive rectangular straight wing, built on a NACA 0012 airfoil. The Reynolds number, made on the upstream flow velocity U ∞ and the airfoil chord c, is 30, 000. A sinusoidal plunging motion along the vertical z axis is considered [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] :

h = h max cos 2π f t ( 1
)
where h is the wing position along the z axis, h max = 0.175c its maximum value, f the plunging motion frequency and t the time. A comparison between the reduced order LDVM, the dynamic stall model [START_REF] Bøckmann | Wave propulsion of ships[END_REF] and the Unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation is presented in 2D, with the development of the propulsive coefficients with respect to the reduced flapping frequency. Then, an extension of the 3D LDVM [START_REF] Faure | High angle-of-attack aerodynamics of a straight wing with finite span using a discrete vortex method[END_REF] is proposed in order to take into account wing dihedral, resulting from the spanwise flexibility. The 3D configuration is also considered using a pseudo-stationary lifting line theory and the 3D LDVM, for the aforementioned three wings with different stiffnesses, referred to as inflexible, flexible and highly flexible cases. The thrust, lift and quarter-chord pitching moment coefficients are given for these flexibility cases.

The comparison between LDVM predictions and experimental data of the vortices shed from the wing is presented.

II. NUMERICAL METHODS

A. Dynamic stall model and pseudo-stationary non-linear lifting line theory

Dynamic stall model

The 2D model used hereafter is a low computational cost model presented by Bøckmann [START_REF] Bøckmann | Wave propulsion of ships[END_REF] , which derives from the Beddoes-Leishman dynamic stall model [START_REF] Leishman | A semi-empirical model for dynamic stall[END_REF] . It provides the unsteady aerodynamic forces acting on a 2D airfoil with an arbitrary development of the angle of attack α. The model is composed of three main contributions representing various flow regimes: unsteady attached flow, trailing edge separation and leading edge separation. Some empirical time lags are included to model the vortex dynamics properly. In particular, the lift due to a vortex developing near the leading edge is also accounted for. Bøckmann adds some modifications: in particular, the effect of compressibility is removed. The non-circulatory normal force (also known as added mass force) is calculated differently and as proposed in Moriarty & Hansen [START_REF] Moriarty | Aerodyn theory manual[END_REF] , the separation point locations are expressed slightly differently, in order to be able to recreate the static data, for large values of the angle of attack. Static forces coefficients tables are used to calculate the separation points location. In present study, in order to be consistent with the experiment [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] , static lift data [START_REF] Cleaver | Vortex mode bifurcation and lift force of a plunging airfoil at low reynolds numbers[END_REF] are used for a Reynolds number of 30, 000 in the range -20 • ≤ α ≤ 20 • . For higher angles of attack and static drag, data [START_REF] Sheldahl | Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines[END_REF] for a Reynolds number of 20, 000 are used. The method used in 3D is a pseudo-stationary version of the non-linear lifting-line theory of Prandtl [START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] . The non-linear lifting-line theory is a reduced order method giving the lift and drag on a finite wing with sweep and/or dihedral in an incompressible, irrotational and steady flow [START_REF] Duport | Benchmarking of a 3D non-linear lifting line method against 3D RANSE simulations[END_REF][START_REF] Duport | Comparison of 3D non-linear lifting line method calculations with 3D RANSE simulations and application to the prediction of the global loading on a cornering kite[END_REF][START_REF] Phillips | Modern adaptation of prandtl's classic lifting-line theory[END_REF] . The wing is sliced into sections with constant circulations. 3D effects are modeled by a vortex sheet composed of horseshoe vortices starting at the quarter chord of the wing (figure 1). Solving the non-linear problem consists in finding the piece-wise constant circulation evolution along the span, satisfying both Kutta-Joukowski theorem and 2D external data for each section. The external data are 2D lift and drag curves relatives to the angle of attack, coming from experiments or higher order simulations, and in that sense viscosity can be accounted for. In present work, a thin symmetrical airfoil law is used, C L = 2πα. The wing oscillation is modeled by a succession of stationary calculations. At each time step, the shape of the wing and the apparent velocities of the sections are updated and the non-linear problem is solved. The aerodynamic forces are then integrated over the oscillation cycle and divided by the period to obtain the mean aerodynamic coefficients.

B. Reynolds-averaged numerical simulation

2D URANS are carried out for a Reynolds number of 30, 000, following a finite volume method using the Open-Foam software suite. Solver parameters, mesh and turbulence modeling are validated on static cases, compared to literature results. The numerical scheme is a merge between Pressure Implicit with Splitting of Operators (PISO) and Semi-Implicit Method for Pressure Linked Equation (SIMPLE), second order in space and implicit Euler (first order) in time. Boundary layer resolution (see below) and the need to predict a moving laminar separation followed by prospective turbulent reattachment, lead to a refined mesh around the airfoil and, through stability considerations, a very short time step (the order of 5 × 10 -5 / f ).

A C-shape structured mesh is used, and independency checked by monitoring the values of aerodynamic coefficients and location of characteristic stations (separation and reattachment for the laminar separation bubble obtained in static computations), for mesh sizes ranging from 82, 370 to 296, 300 points. I) and laminar separation bubble position are concerned. So that the 163, 000 mesh points grid is kept for the dynamic simulations. These values are to be compared with the drag coefficient value of 0.0218 obtained by Galbraith & Visbal [START_REF] Galbraith | Implicit large eddy simulation of low reynolds number flow past the sd7003 airfoil[END_REF] , who used high-order LES to draw the SD 7003 polar plot at the same Reynolds number, and the LES result of 0.0226 obtained by Cataleno & Tognaccini [START_REF] Catalano | Influence of free-stream turbulence on simulations of laminar separation bubbles[END_REF] . During dynamic runs, the mesh is deforming to follow the movement of the airfoil. Deformation is computed using Laplacian diffusion of mesh points displacement and preserves the mesh resolution in the vicinity of the wall. Turbulence modeling is given special attention as the flow is likely to be transitional on the airfoil at the considered Reynolds number. The present work follows the idea, initially exposed by Spalart & Strelets [START_REF] Spalart | Mechanisms of transition and heat transfer in a separation bubble[END_REF] and further exemplified by Crivellini & d'Alessandro 28 , according to which, with some care, the Spalart-Allmaras model can correctly predict the flow around an airfoil at low Reynolds number, including the so-called laminar separation bubble and induced transition to turbulence. As exposed in the aforementioned papers, the mesh is designed so as to ensure that wall ordinate remains the order of unity on the whole wall, and inflow turbulent viscosity ν t is set at a very low level. The value ν t = 10 -14 ν is selected, with ν the flow viscosity. Results published for a SD 7003 airfoil at 4 • of angle of attack and chord Reynolds number of 60, 000 are used for further validation [START_REF] Ol | Comparison of laminar separation bubble measurements on a low Reynolds number airfoil in three facilities[END_REF] . Simulations are carried out for the same airfoil, angle of attack and Reynolds number, using the retained mesh and solver parameters, and are found to correctly predict the location of separation and reattachment of the laminar separation bubble: the two finest meshes lead to the same values of 0.28 chord length and 0.66 chord length for the separation and reattachment respectively, while results compiled in Ol et al. [START_REF] Ol | Comparison of laminar separation bubble measurements on a low Reynolds number airfoil in three facilities[END_REF] ranged from 0.2 to 0.33 for separation and 0.59 to 0.63 for reattachment.

C. LDVM 1. 2D LDVM

The 2D LDVM is a reduced order unsteady flow simulation method around an airfoil which requires low computing power and short simulation time [START_REF] Ramesh | Discrete-vortex method with novel shedding criterion for unsteady aerofoil flows with intermittent leading-edge vortex shedding[END_REF][START_REF] Faure | A modified discrete-vortex method algorithm with shedding criterion for aerodynamic coefficients prediction at high angle of attack[END_REF][START_REF] Faure | Numerical study of two-airfoil arrangements by a discrete vortex method[END_REF] . It is based on the thin airfoil theory and models the time-dependent vorticity distribution along the chord-wise x axis as a Fourier series [START_REF] Glauert | The Elements of Aerofoil and Airscrew Theory[END_REF] :

γ (θ ,t) = 2U ∞ A 0 (t) 1 + cos θ sin θ + ∞ ∑ n=1 A n (t) sin nθ (2)
where U ∞ is the upstream flow velocity, A 0 , . . . , A n are the time-dependent Fourier coefficients calculated from the induced velocity w i normal to the airfoil [START_REF] Katz | Low-Speed Aerodynamics[END_REF] and the new variable θ resulting from the transformation of the chord-wise coordinate x such as:

x = c 2 (1 -cos θ ) (3)
The bound circulation attached to the airfoil is then:

Γ B (t) = ¢ c 0 γ(x,t)dx (4)
The method considers an unsteady potential flow augmented by the modeling of a possible detachment, according to the Leading Edge Separation Parameter (LESP), which reaches a threshold (LESP crit ) if a leading edge separation occurs. The LESP and its critical value can be obtained with a technique, which uses a few pressures in the airfoil leading edge region for deducing the aerodynamic state of an airfoil [START_REF] Saini | Leading-edge flow sensing for detection of vortex shedding from airfoils in unsteady flows[END_REF] . In this unsteady method, for each time step, a Trailing Edge Vortex (TEV) is computed if the flow is attached. If the flow is detached, according to the LESP, in addition to the TEV, a Leading Edge Vortex (LEV) is computed. For further details on the LDVM, the reader can refer to Faure et al. [START_REF] Faure | A modified discrete-vortex method algorithm with shedding criterion for aerodynamic coefficients prediction at high angle of attack[END_REF] . For the NACA 0012 airfoil considered hereafter [START_REF] Chiereghin | Unsteady lift and moment of a periodically plunging airfoil[END_REF] , the value of LESP crit is 0.25 for the Reynolds number range 10, 000 -30, 000. The viscosity is a key feature determining the optimum oscillation frequency of an airfoil [START_REF] Floryan | Efficient cruising for swimming and flying animals is dictated by fluid drag[END_REF] . Thus, in present study, a viscous drag correction is added when the flow is attached, in order to stray away from the purely potential thrust prediction. This situation arises either for k G < 1.2 during the whole cycle, or during part of the cycle for higher frequencies. The correction consists in the addition, to the computed drag coefficient, of the static drag coefficient C D (α) in the attached flow region [START_REF] Sheldahl | Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines[END_REF] . This drag coefficient is obtained considering the relative angle of attack α rel , which is the angle between the chord and the combination of the upstream velocity U ∞ and the relative velocity induced by the heave motion ḣ, for each time step:

α rel = arctan ḣ U ∞ (5)
with ˙the time derivation. If the flow is detached (LESP > LESP crit ), no more correction is added since the drag is mainly due to flow detachment.

3D LDVM

A coupling between the 2D LDVM and the lifting line theory [START_REF] Prandtl | Applications of modern hydrodynamics to aeronautics[END_REF] has been developed for a finite aspect ratio unswept wing [START_REF] Faure | High angle-of-attack aerodynamics of a straight wing with finite span using a discrete vortex method[END_REF][START_REF] Ramesh | Experimental and numerical investigation of post-flutter limit cycle oscillations on a cantilevered flat plate[END_REF] , consisting in 2D LDVM calculations in each spanwise sections, with a coupling with the lifting line theory. This approach results in the addition of a correction term A 0,3D in the formulation of the A 0 Fourier coefficient. Consider a wing of finite span b and chord c with no sweep or dihedral, decomposed into N sections along its span. The span direction is y, with the variable ψ such as:

y = - b 2 cos ψ (6)
For further details on 3D LDVM applied to a straight wing, the reader can refer to Faure & Leogrande [START_REF] Faure | High angle-of-attack aerodynamics of a straight wing with finite span using a discrete vortex method[END_REF] . For a rigid wing with no sweep, the vortices representing the spanwise vorticity distribution are collinear and their mutual interaction is null. But this is no more true for a wing with dihedral, induced for instance, by spanwise flexibility as in present study. Thus, in addition to the 3D LDVM, to take into account a possible dihedral [START_REF] Phillips | Modern adaptation of prandtl's classic lifting-line theory[END_REF] , the inclusion of a weak coupling between the different wing sections is implemented in the method (figure 2). It is based on the interaction between the bound vorticity Γ B, j , Γ B,l of different spanwise sections of the wing [START_REF] Duport | Comparison of 3D non-linear lifting line method calculations with 3D RANSE simulations and application to the prediction of the global loading on a cornering kite[END_REF] , but neglects the interaction between the shed vortices developing in different spanwise sections. However, as in 2D LDVM, there are interactions between LEV, TEV and bound vorticity in a given section. The velocity normal to the airfoil section results in the addition of a coupling term between the other sections:

w i (y,t) = w i,rigid (y,t) + w i, f lex (y,t) (7) = -U ∞ N ∑ n=1 nP n sin nψ sin ψ + 1 4π N ∑ j=1 Γ B, j (ψ) N ∑ l=1,l = j 1 r 2 l, j
where w i,rigid is the induced velocity resulting from the cou-FIG. 2: Sketch of the weak coupling in the 3D LDVM: the interaction between the different spanwise sections is considered for the bound circulation (in green), for a wing with dihedral.

pling between 2D LDVM and lifting line theory for a straight and rigid wing, w i, f lex is the induced velocity resulting from the interaction between the N wing sections due to dihedral, and P n the coefficients of the lifting-line theory such as:

Γ B (ψ,t) = 2bU ∞ N ∑ n=1 P n (t) sin nψ (8)
with:

r l, j = y f ,l -y f , j 2 + z f ,l -z f , j 2 (9)
with (y f ,l , z f ,l ) and (y f , j , z f , j ) the coordinates of wing sections l and j. As a result, the modification of the 3D LDVM for a wing with dihedral appears in the expression of the correction term:

A 0,3D (t) = N ∑ n=1 nP n sin nψ sin ψ - 1 4πU ∞ N ∑ j=1 Γ B, j N ∑ l=1,l = j 1 r 2 l, j (10)
The time-dependent Fourier coefficients of the thin airfoil theory (Eq. 2) are written as [START_REF] Faure | High angle-of-attack aerodynamics of a straight wing with finite span using a discrete vortex method[END_REF][START_REF] Ramesh | Augmentation of inviscid airfoil theory to predict and model 2D unsteady vortex dominated flows[END_REF] :

A 0 (t) = A 0,2D (t) + A 0,3D (t) (11) A n (t) = A n,2D (t)
with A 0,2D , . . . , A n,2D are the time-dependent Fourier coefficients for the 2D LDVM defined in equation 2.

III. 2D RESULTS

The Garrick reduced frequency [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF][START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] is defined as:

k G = π f c U ∞ (12)
In this section relative to the 2D flow, different reduced frequencies between 0.5 and 5.5 are considered. The maximum heave amplitude of the wing is h max = 0.175c. The thrust T h and power P coefficients are defined as:

C T = T h 1 2 ρU 2 ∞ c (13)
C P = L ḣ 1 2 ρU 3 ∞ c ( 14
)
with ρ the fluid density, L the lift and ḣ the vertical wing root velocity. The propulsive efficiency is obtained as:

η = C T C P ( 15
)
where denotes time-averaging over a complete number of oscillations. The thrust coefficient increases with the reduced frequency (figure 3), a good agreement is found between the simulations, dynamic stall model, experiments and LDVM, while much larger values are found for the potential methods which fail to predict the flow detachment. For reduced frequencies lower than 1.2, no flow detachment is observed in the LDVM which behaves in a way similar to the potential method, and for 1.2 ≤ k G ≤ 1.5, the detached region is relatively limited during a cycle. A small bump is observed in the LDVM prediction of C T for k G ≤ 1.7 (figure 3), and is understood as a drawback of the method to correctly predict the viscous drag. That is the reason why a viscous correction has been implemented when the flow is attached. Similar agreement is observed for the power coefficient (figure 4), with a lower estimate for the potential methods and dynamic stall model. Note that, according to the experimental results, there is very little effect of the Reynolds number in the range of study. The propulsive efficiency (figure 5) shows a maximum for reduced frequencies between 1 and 1.3, with maximum values between 0.27 and 0.4. Some discrepancies are observed between the references, because of the different approach and the low numerical or experimental signal-to-noise ratio for low frequencies, with a laminar or transitional flow. As a consequence, for this low frequency range, the efficiency is the ratio of two small quantities, which can explain the dispersion. The dynamic stall model predicts the efficiency, within the range of other methods, up to k G = 2 but gets a threshold around η = 0.28 for larger reduced frequencies, because of the under-prediction of the power coefficient. This may be caused by the calibration of the dynamic stall model which is conducted for k G ≤ 0.5 [START_REF] Bøckmann | Wave propulsion of ships[END_REF] . When the frequency is low k G < 1.5, the results of the LDVM without viscous correction agree with the panel method, whereas with the viscous correction there is a better agreement with the experimental and other numerical results.

The URANS points are slightly overestimating the predictions of previous simulation 9 and experiment [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] . A good agreement is found for k G = 1.82 between dynamic stall model, Navier-Stokes solver [START_REF] Young | Oscillation frequency and amplitude effects on the wake of a plunging airfoil[END_REF] , URANS and LDVM with viscous correction. All these methods present values of efficiency approximately 25% higher than the experimental data [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] .

Figure 6 presents the comparison of the 2D simulations between URANS and LDVM at a Reynolds number of 30, 000 for four different relative time steps over a period T of the plunge. The following flow fields are obtained after at least one flapping cycle. Note the general feature of the flow de-veloping downstream of the airfoil, with a typical reverse Von Kármán vortex street, with an anticlockwise vortex row above the airfoil chord axis and a clockwise vortex row below the airfoil chord axis, typical of a thrust producing wake [START_REF] Young | Oscillation frequency and amplitude effects on the wake of a plunging airfoil[END_REF] . For each relative time, over a plunging period, the positions of the shed vortices are in accordance between the two methods. Most of the vortices generated on the suction and pressure sides of the airfoil, are also equally depicted. Despite the difference between the simulation method, no noticeable change in the flow feature is observed. For t/T = 0 (figures 6a and 6e), the airfoil reaches its maximum amplitude and its velocity is null. Negative vorticity vortices are advected on the suction side towards the trailing edge while a positive vorticity vortex is present on the pressure side and a vortex is starting shedding from the trailing edge. For t/T = 0.25 (figures 6b and 6f), the airfoil gets its maximum downward vertical velocity. The pressure side positive vorticity vortex is advected and present around the mid-chord, the negative vorticity vortices are just before the trailing edge and the positive vortex develops downstream from the airfoil. A new cluster of negative vorticity is starting near the suction side leading edge. For t/T = 0.5 (figures 6c and 6g), the airfoil reaches its minimum amplitude and its velocity is null again. Negative vorticity near the leading edge develops on the suction side, a negative vorticity vortex is shedding from the trailing edge and the pressure side positive vorticity is advected and placed after the mid-chord. For t/T = 0.75 (figures 6d and 6h), the airfoil achieves its maximum upward vertical velocity. The shed negative vortex is developing while the previously shed positive vortex is completely detached.

IV. 3D RESULTS

A. Spanwise wing deformation

The fluid density and spanwise flexibility play an important role in the steady and dynamic performance of the wing [START_REF] Chae | Influence of spanwise flexibility on steady and dynamic responses of airfoils vs hydrofoils[END_REF] . For a flexible wing of span b, the vertical deformation z as a function of the span direction y and time t, is obtained by the formulation [START_REF] Liu | Propulsive performance from oscillating propulsors with spanwise flexibility[END_REF] :

z (y,t) = h (t) + d (y,t) (16)
with h (t) the wing root sinusoidal plunging motion along the vertical axis z and d (y,t) the function governing the shape of the spanwise deflection:

d (y,t) = d max 2 2y b 2 ∓ 4 3 2y b 3 + 1 3 2y b 4 sin (ωt + Φ) (17)
with d max the maximum deflection at the wing tip, ω the heave angular frequency and Φ the spanwise deflection phase angle relative to the heave motion. In present formulation, the maximum deformation is located at the wing tip. The values of the parameters d 0 and Φ are determined from the inflexible, flexible and highly flexible experimental cases [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] and converted into the formulation [START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] [START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] for the three wings of different stiffnesses and k G = 1.82.

The relative vertical position of the wing root and tip z/c for the oscillatory motion in pure heave at a Reynolds number of 30, 000 is given in figure 7 for two periods [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] . The reduced frequency is k G = 1.82. For the 3D study, the wing aspect ratio considered is A R = 6, taking into account the side root wall mirror effect. The continuous black line is the vertical displacement of the wing root. For the rigid wing, the same motion is theoretically expected at the wing tip. However, because of some limited flexibility in the wing structure, the observed vertical motion of the tip in the experiment is plotted by the dotted green line, and the case is hereafter referred to as the inflexible case. For the flexible and highly flexible wings, larger tip amplitudes and phase lags are observed.

B. Inexible wing

In this section, we compare experimental, numerical and theoretical results [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF][START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF][START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] with 2D LDVM and 3D LDVM simulations, for an inflexible wing. The thrust coefficient is presented for two periods of the vertical oscillation in figure 8. A two peak pattern is observed for each cycle, which is consistent with the shedding of two vortices downstream from the airfoil. A general agreement is found between the available data and the 2D or 3D LDVM. The 2D estimates (URANS and LDVM 2D) slightly overestimate the 3D predictions (Navier-Stokes solver of Gordnier et al. [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] and LDVM 3D), while the measurements of Heathcote et al. [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] present a lower peak of negative thrust for each period around t/T = 0.55, which may be due to the experimental device. The approach based on the pseudo-stationary non-linear lifting line theory [START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] provides larger maximum peak values.

Figure 9 shows the lift coefficient. This quantity is one order of magnitude larger than the thrust coefficient, and a good agreement is found between all the 2D or 3D data sets. Similar comments are valid for the quarter chord pitching moment coefficient (figure 10) with some discrepancies for 2D LDVM after the first cycle. The larger amplitude observed in 2D LDVM after the first period could be caused by the influence, on the pitching moment, of a larger shed vortex created after the first period. This vortex could be tempered in 3D LDVM by the spanwise effects imposed by the lifting line correction.

C. Flexible wing

For the flexible wing, a similar two peak pattern is observed for each cycle in the thrust coefficient (figure 11). Despite some dispersion between the maximum peak values, note that in each data set, the thrust coefficient is larger than what is found for the inflexible wing, indicating a beneficial effect of a limited flexibility. For the lift and moment coefficients (figures 12 and 13), there is few published data and only limited possible comparison [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] showing a relatively good agreement, with the pseudo-stationary lifting line theory (for the lift), and 3D LDVM.

D. Highly exible wing

The thrust coefficient for the highly flexible wing is given in figure 14. As previously observed [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF][START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF][START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] , the thrust coeffi- cient is much lower than for the inflexible or flexible wing. A two peak per cycle pattern is observed in some data [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF][START_REF] Simonet | Étude de l'efficience d'un système de propulsion maritime par foil oscillant souple[END_REF] but there is no clear periodic pattern in measurements of Heathcote et al. [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] and 3D LDVM. For the lift and moment coefficients (figures 15 and 16) the amplitudes are comparable but a surprising phase shift of half a period is observed between the high-fidelity aeroelastic simulation of Gornier et al. [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] and the pseudo-stationary lifting line theory model and 3D LDVM. It may be due to a second order flexion mode of the highly flexible wing, resulting from the results of the structural dynamics solver [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] , which is not considered in the formulation of the spanwise deflection (equation 17). (PIV) measurements with LDVM results. The upstream flow is in the positive direction of the x axis, and to minimize the overlap of the fields, the y axis is scaled differently to the x and z axes. The six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the PIV measurements and LDVM simulations. In these figures, the arrows indicate the direction of motion of the wing root and tip. For PIV measurements, only the flow downstream from the trailing edge is plotted, while the suction and pressure side flows are also plotted for LDVM. In addition, in LDVM flow fields, the rigid wing position is represented by light dotted lines and the actual (flexible) wing position by bold dotted lines. For t/T = 0, in the first section from the root (y/c = 0.5), two vortices of opposite sign are visible (figures 17a and 17b), and the flow is similar to the 2D flow (figure 6). In this plane, the flow pattern is identical for the three wings with a stronger vorticity for the flexible wing and a weaker intensity for the highly flexible wing. The vorticity pattern is similar in section y/c = 1, but nearer the tip for y/c ≥ 2, the clockwise vortex (light shading or red points) forks into two clusters. In LDVM, this behavior is observed except for the inflexible wing. For the highly flexible wing (figure 18b), the vorticity centers for y/c ≥ 1.5, which are forming larger scattered clusters, are the evidence of a more complex shedding downstream from the trailing edge.

E. Comparison between the three wings

For t/T = 0.25 (figure 18), corresponding to a downward motion of the wing root, the tips of the inflexible and flexible wings move in the same direction, but the tip of the highly flexible wing moves in the opposite direction. This results in vorticity of opposite sign shed in the flow, immediately downstream from the trailing edge, between the tip and the root, for the highly flexible wing. Thus, positive vorticity (light shade or red points) is created immediately downstream from the trailing edge for y/c ≤ 2, while negative vorticity (black shade or blue points) is created immediately downstream from the trailing edge for y/c > 2. The wing flexibility is an important parameter for the shed vortices, and opposite vorticity generation between wing root and tip is correlated with lower [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] with LDVM simulation for the very same stiffness parameters presented in table II and a similar point of view of the 3D wing, with the wing root in the foreground and the tip in the background. Note that the six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the computations downstream of the trailing edge [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] and LDVM simulations. Because of a longer time simulated and a larger downstream view, an additional downstream vortex is plotted in LDVM. This vortex, located farther from the trailing edge, is not depicted in the high fidelity computations of the vorticity field. For the flexible case and t/T = 0, the similarity between high fidelity computation and LDVM is observed in the vortex shedding from the trailing edge (figures 19a and 19c) with a rather 2D flow up to y/c = 2. The same comment is valid for t/T = 0.25, with the division of the negative vorticity structure into several scattered vortices near the wing tip (figures 20a and 20c). For the highly flexible case, it has been reported that the differences in comparison between experiments [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] and computations [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] are attributed in part to the phase shift between the structural deflection. This discrepancy is probably due to the deflection formulation (Eq. 17) which does not take into account the higher order deflection modes obtained by the structural solver. However, if the computed deflection parameters obtained by Gordnier et al. [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] and presented in table II are used in LDVM simulation, a good agreement is found between the shed vortices (for t/T = 0 in figures 19b and 19d and for t/T = 0.25 in figures 20b and 20d).

CONCLUSION

A low order discrete vortex method is considered to predict the performance and flow field produced by a periodic flapping wing. The LDVM is validated on a benchmark configuration of the sinusoidal heave motion of a rectangular wing of aspect ratio of 6, for Reynolds numbers of the order of 30, 000. First, a 2D characterization of the thrust and power coefficients, propulsive efficiency and flow field is conducted, comparing low order dynamic stall model and LDVM with a URANS simulation. A good agreement is found between the last two methods and published data. Then, 3D pseudostationary lifting line theory and LDVM are adapted to consider the wing spanwise flexibility. The integration of a dihedral, resulting from the spanwise wing deformation, is integrated in the method. Three wings with different stiffnesses are considered from previous experiment [START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] : inflexible, moderately flexible and highly flexible cases. A general agreement between the methods is found for the thrust, lift and quarter chord pitching moment. Note that, as previously published, the intermediate flexible wing provides a better thrust coefficient while excessive flexibility proves to be detrimental. For the highly flexible wing, there is a larger dispersion between simulations and experiment for lift coefficient. In particular, the lift and moment coefficients show a phase shift of half a period between the numerical simulation [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] and the pseudostationary lifting line theory and 3D LDVM. This may be due to higher order flexion modes present for the highly flexible wing. These modes could result from the structural dynamics solver [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] , which is not considered in the present formulation of the spanwise deflection. For the inflexible and flexible cases, comparisons between measured flow features and LDVM predictions show the development of similar vortex clusters during the oscillating motion. The position of the shed vortices from the vorticity field is accurately predicted by the method for the highly flexible case, whenever proper stiffness is considered [START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] . The main result of the study is that a low order 3D discrete vortex method can predict the aerodynamic coefficients and the flow features of the periodic heave motion of wings of different stiffnesses, with a good accuracy and at a small computing cost. The LDVM is able to run on a laptop. Reported efficiencies between 83% and 92% are found [START_REF] Liu | Propulsive performance from oscillating propulsors with spanwise flexibility[END_REF][START_REF] Fish | Comparative kinematics and hydrodynamics of odontocete cetaceans: morphological and ecological correlates with swimming performance[END_REF][START_REF] Fish | Hydrodynamic performance of aquatic flapping: efficiency of underwater flight in the manta[END_REF][START_REF] Ayancik | Three-dimensional scaling laws of cetacean propulsion characterize the hydrodynamic interplay of flukes' shape and kinematics[END_REF] for the propulsion by periodic flapping of a fluke used by the cetaceans. These wing flows are developing for much larger Reynolds numbers and amplitudes, combining heave and pitch, maintaining the wing section relative angle of attack to a small value, with an attached flow. This configuration could be investigated in a similar way from present study using the LDVM.

FIG. 1 :

 1 FIG. 1: Example of a low discretized lifting line model with sweep and dihedral 22 with resulting lift L and thrust T h .

Figures 3 to 5 FIG. 4 :

 54 Figures 3 to 5 present the comparison of the thrust and power coefficients and the propulsive efficiency for different 2D references: analytical model 6 , numerical panel method, flow simulation using a compressible Navier-Stokes solver[START_REF] Young | Oscillation frequency and amplitude effects on the wake of a plunging airfoil[END_REF][START_REF] Young | Numerical simulation of the unsteady aerodynamics of flapping airfoils[END_REF] , water tunnel experiments[START_REF] Heathcote | Effect of spanwise flexibility on flapping wing propulsion[END_REF] and present dynamic stall model, URANS and LDVM.

FIG. 5 :

 5 FIG. 5: Comparison of propulsive efficiency versus Garrick frequency k G for a NACA0012 airfoil in pure heave, for present study dynamic stall model, LDVM and URANS with published data of linear theory (Garrick 6 ), panel method, Navier-Stokes solver (Young & Lai 9 ) and experiment (Heathcote et al. 11 ).

FIG. 6 :FIG. 7 :

 67 FIG. 6: NACA 0012 airfoil in pure heave for a Reynolds number of 30, 000: comparison between 2D finite volume simulation (first row) and 2D LDVM simulation (second row) for different relative times. The dimensionless vorticity Ωc/U ∞ is provided for finite volume simulation: (a) t/T = 0, (b) t/T = 0.25, (c) t/T = 0.5, (d) t/T = 0.75 and the dimensionless circulation Γ/(U ∞ c) is represented for LDVM : (e) t/T = 0, ( f ) t/T = 0.25, (g) t/T = 0.5, (h) t/T = 0.75.

FIG. 8 :

 8 FIG. 8: Comparison the inflexible wing, of the time histories of the thrust coefficient from experimental results of Heathcote et al. 11 and computational results of Gordnier et al. 12 with present study results of the pseudo-stationary lifting line theory, 2D URANS, 2D LDVM and 3D LDVM.

FIG. 9 :

 9 FIG.9: Comparison for the inflexible wing, of the time histories of the lift coefficient from computational results of Gordnier et al.[START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] with present study results of the pseudo-stationary lifting line theory, 2D URANS, 2D LDVM and 3D LDVM.

Figures 17 and 18

 18 Figures 17 and 18 represent the fields of the y axis vorticity component Ω at two relative times for the three wings of different stiffnesses, and compare Particle Image Velocimetry

FIG. 13 :

 13 FIG.13: Comparison for the flexible wing, of the time histories of the quarter-chord pitching moment coefficient from computational results of Gordnier et al.[START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] with present study results of the 3D LDVM.

FIG. 15 :

 15 FIG.15: Comparison for the highly flexible wing, of the time histories of the lift coefficient from computational results of Gordnier et al.[START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] with present study results of the pseudo-stationary lifting line theory and 3D LDVM.

FIG. 17 :

 17 FIG. 17: Spanwise vorticity fields for t/T = 0 for the inflexible (top), flexible (middle), and highly flexible (bottom) wings; Re = 30, 000, k G = 1.82: (a) dimensionless vorticity Ωc/U ∞ contours plotted from PIV measurements downstream of the trailing edge 11 Reproduced with permission from J. Fluids Struct. 24 (2008). Copyright 2008 Elsevier, (b) dimensionless circulation Γ/(U ∞ c) vortical centers from LDVM simulation. The six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the PIV measurements and LDVM simulations.

FIG. 18 :

 18 FIG. 18: Spanwise vorticity fields for t/T = 0.25 for the inflexible (top), flexible (middle), and highly flexible (bottom) wings; Re = 30, 000, k G = 1.82: (a) dimensionless vorticity Ωc/U ∞ contours plotted from PIV measurements downstream of the trailing edge 11 Reproduced with permission from J. Fluids Struct. 24 (2008). Copyright 2008 Elsevier, (b) dimensionless circulation Γ/(U ∞ c) vortical centers from LDVM simulation. The six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the PIV measurements and LDVM simulations.

FIG. 19 :

 19 FIG. 19: Spanwise vorticity fields for t/T = 0 for the flexible (top) and highly flexible (bottom) wings; Re = 30, 000, k G = 1.82: dimensionless vorticity Ωc/U ∞ contours plotted from computations downstream of the trailing edge 12 : (a) flexible, Reproduced with permission from J. Fluids Struct. 40 (2013). Copyright 2013 Elsevier, (b) highly flexible, Reproduced with permission from J. Fluids Struct. 40 (2013). Copyright 2013 Elsevier, dimensionless circulation Γ/(U ∞ c) vortical centers from LDVM simulation: (c) flexible, (d) highly flexible. The six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the computations downstream of the trailing edge 12 and LDVM simulations.

FIG. 20 :

 20 FIG. 20: Spanwise vorticity fields for t/T = 0.25 for the flexible (top) and highly flexible (bottom) wings; Re = 30, 000, k G = 1.82: dimensionless vorticity Ωc/U ∞ contours plotted from computations downstream of the trailing edge 12 : (a) flexible, Reproduced with permission from J. Fluids Struct. 40 (2013). Copyright 2013 Elsevier, (b) highly flexible, Reproduced with permission from J. Fluids Struct. 40 (2013). Copyright 2013 Elsevier, dimensionless circulation Γ/(U ∞ c) vortical centers from LDVM simulation: (c) flexible, (d) highly flexible. The six spanwise visualization sections (y/c =0.5, 1, 1.5, 2, 2.5, 2.85) are identical for the computations downstream of the trailing edge 12 and LDVM simulations.

TABLE I

 I

	: Lift, drag and quarter chord pitching moment
	coefficients for an angle of attack of 4 • and different meshes.
	296, 300, are found to give the same results as far as global
	coefficients values (table

  of equation 17 and presented in table II. Note that the computation of Gordnier et al.[START_REF] Gordnier | Highfidelity aeroelastic computations of a flapping wing with spanwise flexibility[END_REF] includes a structural dynamic solver, and some phase shift is observed between the experimental values, especially for the highly flexible case.

	wing stiffness	d 0 /c	Φ
	inflexible Heathcote et al. (2008) 0.0265 56.9 •
	flexible	Heathcote et al. (2008) 0.148 34.9 •
		Gordnier et al. (2013) 0.146 29.68 •
	highly flexible Heathcote et al. (2008) 0.418 -48.9 •
		Gordnier et al. (2013) 0.444 -61.2 •

TABLE II :

 II Deformation parameters

ACKNOWLEDGEMENTS

The authors are gratefully to Prof. I. Gursul and Dr. Z. Wang for the use of the flow plots in figures 17 and 18 and to Prof. C. S. Cesnik and Prof. P. J. Attar for the use of the vorticity plots in figures 19 and 20.

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article.