Alessandro Fiumara
email: afiumara@assystem.com

Nicolas Gourdain

Vincent Chapin

Julien Senter

Yannick Bury

Numerical and experimental analysis of the flow around a two-element wingsail at Reynolds number 0.53 × 10 6

Keywords:

The rigid wingsail is a propulsion system, utilized in sailing competitions in order to enhance the yacht performance in both upwind and downwind conditions. Nevertheless, this new rig is sensitive to upstream flow variations, making its steering difficult. This issue suggests the need to perform a study on wingsail aerodynamics. Thus this paper reports some investigations done to better understand the flow physics around a scaled model of an America's Cup wingsail, based on a two-element AC72 profile. First a wind tunnel test campaign was carried out to generate a database for aerodynamic phenomena analyses and CFD validation. Unsteady RANS simulations were performed to predict and validate the flow characteristics on the wingsail, in the wind tunnel test conditions. The wind tunnel domain was fully modeled, in order to take into account the facility confinement effects. Numerical simulations in freestream and wind tunnel conditions were then compared with experimental data. This analysis shows the necessity to consider the wind tunnel walls when experimental and numerical data are compared. Numerical simulations correctly reproduce the flow field for low-to-moderate flow angles. However, discrepancies on the pressure distribution increase when the boundary layer starts to separate from the wingsail. In this regard, the flow generated by the slot between both elements of the wingsail is of paramount importance. This slot flow is analyzed in details through PIV measurements and numerical simulations. While the numerical simulation correctly predicts the jet flow itself, it only partially reproduces the interaction between the jet flow and the main flow, especially at high angle of attacks. More precisely, the numerical simulation fails to predict the correct jet flow trajectory, which affects the lift capabilities of the entire wing. The influence of the wingsail deformation during experimental campaigns has been investigated to explain this behavior.

Introduction

Wingsails are increasingly used in sailing competition to substitute conventional soft sails. This new rig, joined to foils, allows the yacht to achieve better performance. However, at the same time sailors may have problems correctly setting and maneuvering the wingsail in all sailing conditions. Some spectacular and dangerous capsizes occur during the last edition of the America's Cup competition, due to this issue. To date, the global performance envelope of wingsails is not completely understood since the aerodynamic phenomena have not been fully investigated. Moreover, the naval environment introduces some perturbations (like atmospheric boundary layer and upstream turbulence) that should be taken into consideration. Furthermore, the aerodynamic characterization has become increasingly important after the introduction of foils, allowing the catamaran to "fly" on the sea ( Fig. 1 ). Without sufficient hydrodynamic lift, significant yacht deceleration can occur, thereby making the research of stable navigation conditions essential.

Typically, wingsails are composed of a main element and a flap, set in a way to obtain maximum performances on the water. Due to high flap deflection angle variations (from 15 °to 40 °) during the navigation, the flow around a wingsail presents some similarities with the flow around an aeronautical wing in high-lift condition. This analogy explains why sailors have already drawn on the aeronautical know-how to enhance the wingsail performance, like the slotted flap. Unfortunately, the wingsail design imposes some restrictions; for example, to reduce the weight of the wingsail, the flap mechanism is based on a unique rotation while, on aircraft, the flap has more complex kinematics. This constraint reduces the Axes of the wingsail reference system X LSB

x-coordinate of the laminar separation bubble x rot

x-coordinate of the flap rotation axis x v , y v , z v Axes of the wind tunnel reference system y F y position of flap L.E. y + Dimensionless wall distance z * Normalized height position z/H ability to correctly set the slot size. Furthermore, the need to tack from both catamaran sides constrains the wingsail to the usage of only symmetric airfoils which have a lower performance than the asymmetric ones normally used in aeronautics.

Only few experimental works exist today ( e.g. [START_REF] Magherini | Parameters affecting the performance of the C-class wingsail[END_REF][START_REF] Blakeley | Design and optimisation of multi-element wing sails for multihull yachts[END_REF] and a description of flow phenomena is rarely proposed. Blakeley's work has also shown the influence of the flap deflection angle and the slot size on multi-element airfoil performances [START_REF] Blakeley | Evaluation of multi-element wingsail aerodynamics form two-dimensional wind tunnel investigations[END_REF]. Additionally, an exhaustive wind tunnel campaign was performed by [START_REF] Viola | Upwind Sail Aerodynamics: a Pressure Distribution Database for the Validation of Numerical Codes[END_REF] for the aerodynamics characterization of soft sails in upwind conditions, but a comparable analysis on wingsails still does not exist.

To close this gap, an experimental campaign was set on a scale model of an America's Cup AC 72 wingsail in the ISAE-Supaero wind tunnel facility. Oil surface flow visualizations and particle image velocimetry tests were performed during the wind tunnel campaign to describe the characteristics of the boundary layer transition, both on the main and on the flap elements, and to investigate the physics of the flow in the slot. The aim was to provide a description of the flow features over a two-element wingsail, tracking the most sensitive and critical zones in the flowfield and to understand the abilities of the RANS approach to make numer- ical predictions. Both low and high flap deflection angle configurations were analyzed, corresponding respectively to upwind and downwind settings. Reynolds number on the scale model in wind tunnel conditions is 0.53 × 10 6 , 20 times smaller compared to actual AC 72 wingsails during navigation [START_REF] Collie | Aerodynamic design development of AC 72 wings[END_REF]. Nevertheless smaller wingsails are currently used on C-class catamarans ( Re = 0.8 × 10 6 ), in the "little America's Cup" competition.

In the first section of this paper a new methodology is proposed to reproduce the wind tunnel domain and its validation is presented. The second section compares the results of numerical simulations (based on unsteady RANS) in wind tunnel and freestream conditions, with the experimental database. Numerical predictions are then further investigated by comparison with oil flow visualizations, to emphasize the role of laminar to turbulent transition and boundary layer separations. Finally, a physical analysis of the jet flow is done by comparing the numerical velocity scalar maps and the numerical solution with the PIV data and discussed in the third section of the paper.

Experimental methodology

A wingsail scale model of the America's Cup class AC 72 was designed and used for the tests ( Fig. 2 ). It is composed of two elements, the main element and the flap, divided by a slot through which the air can flow. The flap can be set at different angles, pivoting on its axis located at 95% of the main axis. The two ele- ments are composed by NACA symmetrical airfoils, allowing wingsail tacking from both sides.

Pressure ports have been set on three sections of the main element (respectively) located at 25%, 50% and 75% of the wingspan. The pressure sensor used for the measurement is a temperature corrected scan with a + / -5 kPa range and an accuracy of + / -0.15%.

The wind tunnel used for the experimental campaign is the S4 facility owned by "Institut Supérieur de l' Aéronautique et de l' Espace" ISAE-Supaero in Toulouse, an open return wind tunnel with open test section ( Fig. 3 ). The duct has an elliptical shape of 3 m × 2 m. The flow is created by the aspiration of three fan drives of 90 kW each, located at the end of the diffuser (st.5 in Fig. 3 ).

The maximum speed in the duct is 42 m/s. To eliminate low frequency oscillations inside the duct (inherent to such open loop configurations), a gap was created in the first section of the diffuser (st. 4). In doing so, the oscillations are dumped by the creation of a secondary flow, exterior to the diffuser, which recirculates air from the gap to the intake of the diffuser itself (st. 3).

The wingsail model was mounted vertically in the duct (st.2) on a rotating plate that allows adjusting the angle of attack. To reduce the interactions between balance and the wingsail, a disk platform was posed at the base of the wing scale model.

The aerodynamic forces are estimated with a six-component balance. The maximum loads bearable by the balance are 2.40 kN for the longitudinal force, 3.00 kN for the transversal force and 0.50 kNm for the heeling and the pitching moments. Fig. 3 details the two reference systems used in this paper. The first one is the duct system ( x v ,y v ,z v ), whose origin is located at the end of the convergent section in correspondence to the symmetry plane of the duct and at the bottom of the convergent. The x-axis is in the convergent-diffuser direction while the z-axis is directed upwards. The wingsail reference system ( x,y,z ) is translated to the previous one in a way that the origin is located on the leading edge of the wing root section, keeping its position in the symmetry plane of the duct ( x = x V + 2.4 c, y = y v , z = z V + 0.022 H ).

Numerical methodology

The extent of separated regions on wingsail surface may be relevant to prefer Large Eddy Simulations (LES) that are supposed to be the most adapted modeling in these flow conditions. Nevertheless because of the huge computing cost of LES and of the large number of configurations to be studied, it was decided to use unsteady RANS for all numerical simulations. The challenge will be to obtain the best results possible with this model.

Numerical simulations were performed to reproduce the wind tunnel test conditions on the wingsail scaled model. The wingsail geometry was numerically reproduced in both the low and high flap deflection angle configurations recreating at the same time the interface disk on the wing root. Because of the low Reynolds number, based on the mean chord of the wingsail ( Re = 0.53 × 10 5 ), the transition effects have also been considered by the use of the transition model γ -Re ϑ. This model was proposed by Menter et al. (2004 a, b ), based on two transport equations modeling the intermittency factor γ and Re ϑ in turbulent k ω-SST model.

Two approaches were tested for the wingsail environment: 1) a classical approach with freestream conditions (the perturbations induced by the wind tunnel walls are neglected) and 2) an environment integrated approach where the full wind tunnel domain is taken into account.

Wingsail in freestream

The freestream domain considered is a box: length L = 31 c (12 c upstream the wing and 28 c downstream), width l = 40 c and height h = 2 H ( Fig. 4 ), with c the wingsail root chord and H the wingspan. The wingsail is assumed to contact the bottom surface of the box (no root leakage flow).

The reference system has the same characteristics as the one described for the wind tunnel domain, i.e. the origin lies on the leading edge of the main root section, the x-axis has the leadingto-trailing edge direction, the z-axis is directed upward. The entire domain was meshed using a polyhedral mesh with prism layers on the wing surface. The boundary conditions imposed on the box are:

• Velocity inlet: on the inlet, leeward and windward and top surfaces; • Pressure outlet: on the outlet surface;

• Slip wall: on the bottom surface.

The settings chosen for the wingsail are ( Fig. 5 ): SET1 : flap angle δ = 15 °and inlet flow angle α = 0 °; SET2 : flap angle δ = 25 °and inlet flow angle α = 0 °.

The simulations were run using the unsteady RANS approach with the k-ω SST turbulence model and γ -Re ϑ transition model.

Wind tunnel modeling

The boundary condition interaction with the wind tunnel walls for a high-lift configuration is a well known problem, as reported in [START_REF] Rogers | Validation of computed high lift flows with significant wind tunnel effects[END_REF] and [START_REF] Nayani | Numerical study of the high-speed leg of a wind tunnel[END_REF]. Since it is difficult to determine the flow characteristics inside a wind tunnel, directly with measurements, a CFD-based database and experimental database are becoming increasingly necessary in order to have a simple estimation of the aerodynamics and to have an accurate flow description during the wind tunnel tests. This field of investigation is paramount for the fluid dynamics community; recent studies have reported investigations done to correct wind tunnel measurements and improved extrapolation techniques to free flight conditions [START_REF] Melber-Wilkending | Project ForMEx -a new CFD approach for transposition of wind tunnel data towards flight conditions[END_REF][START_REF] Melber-Wilkending | Application of advanced CFD tools for high Reynolds number testing[END_REF][START_REF] Ciobaca | Computational and experimental results in the open test section of the aeroacoustic windtunnel braunschweig[END_REF]. These works showed the wind tunnel influence on the wing, by increasing the effective angle of attack compared to freestream conditions. In the case of soft sail configurations, Viola et al. ( 2013) also suggested taking into account for wind tunnels effects. In this case, a first simulation was carried out to estimate the blockage of the wind tunnel equipped when with the mockup. The extracted velocity field was then used as a boundary condition (avoiding any considerations of full wind tunnel during the numerical simulation of the sail). In the analysis reported in this paper, the numerical modeling of the S4 facility was complicated not only by the blockage of the wind tunnel but also by its elliptical shape and the open wind test section, introducing significant modifications in the flowfield. For this reason, the entire wind tunnel geometry had to be reproduced, as already shown by [START_REF] Fiumara | Comparison of wind tunnel and freestream conditions on the numerical predictions of a flow in a two element wingsail[END_REF] .

The wind tunnel domain was created reproducing at first the convergent and diffuser geometries ( Fig. 3 ). The difficulty is to properly close the duct zone in order to make the domain available for the numerical simulation. Initially the duct was closed using a loft surface from the convergent end section (st. 1 in Fig. 3 ) to the diffuser intake section (st. 3 in Fig. 3 ). However this technique did not lead to a proper reproduction of the wind tunnel flow. An overestimation of the pressure gradients was observed in the rear part of the duct and the jet flow of the duct had the tendency to contract. To overcome these problems, the entire test room was created (framed area in Fig. 3 ) reproducing also the gap existing in the diffuser (st. 4 in Fig. 3 ) to take into account for the external recovery flow. The empty domain was meshed with polyhedral cells with prism layers on the convergent walls only ( Fig. 6 ). The entire mesh is made of 1.3 million cells.

A RANS simulation was then run with STAR-CCM + 9.02 using the turbulence model k-ω SST to model turbulence.

A non-slip condition was used on the convergent surface only in order to account for the effects of the boundary layer on the flow inside the duct. On the remaining surfaces, a slip condition was imposed.

To reproduce the aspiration of the fan drives, a pressure outlet condition was chosen for both the intake of the convergent (st. 0 in Fig. 3 ) and the exit of the diffuser (st. 5). The difference of pressure between the inlet and the outlet were set to modify the pressure value at the outlet in a way to obtain a flow velocity of 20 m/s in the two points where the pressure probes of the real wind tunnel are located. The analysis was run on an Intel I7 processor 2.70 GHz, 32GB of RAM. The time of convergence was about 2 h on 4 cores. This RANS simulation was performed using an empty wind tunnel to compare with experimental data.

Validation of empty wind tunnel simulations

The numerical velocity magnitude (V) distribution of the duct was extracted and compared to the experimental data ( Fig. 7 ). The velocity was normalized with respect to the value assumed to correspond to station 2, at the center of the elliptical section ( i.e. Along the x-axis, the CFD and the experimental data agree well. The difference between the two curves is less than 3% at stations 0 and 3. The tendency of the real S4 is to keep a constant velocity in the first half of the duct and then reduce it significantly in the last quarter of the duct. This loss in velocity is caused by the blockage effect due to the presence of the diffuser intake. In the numerical results, this effect is not reproduced and the distribution is quite constant on the rear part on the duct. In the forward part of the duct, the CFD underestimates the velocity distribution. In the numerical simulation the flow has a favorable pressure gradient for the first part of the duct and zero-pressure gradient in the rear part.

y v = 0, z v = 0.56H).
At st. 2 (where the model was placed), experimental and numerical data are in good agreement along the z direction. Here, the numerical solution matches completely with the experimental results up to z v /H = 0.9. Moving upwards, the CFD solution shows the tendency of the flow to contract, reducing the local velocities on the border of the duct. This loss is estimated to be 3% at z v /H = 1.04, a zone near the wingsail tip, but not a zone that influences the wing directly.

Overall, the analysis of numerical results demonstrates that the numerical wind tunnel is able to reproduce the flow at the mockup location in the real duct.

Wingsail in wind tunnel

Numerical simulations of the wingsail in the wind tunnel were then carried out. The wingsail was considered in the two configurations already tested in the freestream case (SET1 and SET2) and placed at station 2 as in the real case.

The entire domain was meshed using Star-CCM + 9.02 with polyhedra. Prism layers were added on the wingsail, on the disk and on the wind tunnel convergent surface. The layers were set in a way to achieve a normalized distance to the wall y + below 0.5 ( Fig. 8 ) on the wingsail and below 20 on the convergent surface. The choice of a low wall y + on the wingsail surface was derived from the validation tests that had shown the sensitivity of the γ -Re ϑ model to the near wall discretization. In the validation tests for the transition model developed by [START_REF] Suluksna | Correlations for modeling transitional boundary layers under influences of freestream turbulence and pressure gradient[END_REF] [START_REF] Malan | Calibrating the γ -Re ϑ transition model for commercial CFD[END_REF] refers, for the high lift case, to a y + at wall ranging from 0.1 to 0.8. This alternative formulation of the γ -Re ϑ model was proposed by [START_REF] Suluksna | Correlations for modeling transitional boundary layers under influences of freestream turbulence and pressure gradient[END_REF] in order to make up for the lack of the original transport equation of intermittency formulated by Menter et al. (2004a , b) .

) (im- plemented in STAR-CCM + ),
The mesh was refined particularly for the gap between the two elements of the wingsail and the wake region. Refinement was imposed also on the shear layers of the border of the duct. A mesh sensitivity study was performed for the SET2 case. The coarsest mesh (Mesh0) counts 18 Million cells. The mesh was then refined modifying the cell size on the wingsail surface and the polyhedron size in the refined zones. The refinement ratio is respectively 0.85, 0.71 and 0.59 for Mesh1, Mesh2 and Mesh 3 with respect to Mesh0. Simulations were performed in order to extract the lift and drag coefficients of the wingsail at the different grid refinements ( Table 1 ). The lift and the drag coefficients start to converge with Mesh2 having a difference of less than 3% on C L and 1% in C D compared to Mesh3. The final mesh retained for the simulation was Mesh2 (32 million cells).

RANS simulations were run with a k-ω SST turbulence model and activating the γ -Re ϑ. The boundary conditions imposed at the inlet and at the outlet of the wind tunnel are the same used for the empty wind tunnel simulation in a way to keep a flow velocity of 20 m/s in the duct.

Simulations were computed on bi-XeonbE5-2670 Octo processors, 2.60 GHz, 64GB RAM. The computation time was about 6 days on 16 cores. A first convergence was obtained on the aerodynamic coefficient after 40 0 0 iterations. At the same time the pressure distribution over the wingsail and particularly the transition and the laminar bubble zones still presented strong oscillations caused by the unsteady characteristics of the transition phenomena. For this reason the RANS simulations were completed using an unsteady RANS approach, for a total time of 1 s using a time step of 2 × 10 -3 s, corresponding to 50 through flow times (the time needed for a particle to move from the leading edge to the trailing edge).

Results

Comparison of the wingsail results in freestream and in the wind tunnel domain

Wind tunnel numerical simulations (WTNS) and freestream numerical simulations (FSNS) were compared to wind tunnel tests (WTT). The lift and drag aerodynamic coefficients were calculated as well as the C P distribution on the three wingsail reference sections ( i.e. z * = 0.25, 0.50 and 0.75). The error introduced by the balance during experimental tests has to be considered. The tool used during the experimental tests can bear, in fact, loads up to 2.40 kN for the drag and 3.00 kN for the lift. Compared to aerodynamic forces produced by the wingsail in the SET2 configuration ( i.e. D = 0.027 kN, L = 0.20 kN), it represents respectively 1% and 6% of full scale. Therefore, tests on the balance were expressly carried out to give an estimation of the measurement error due to the balance sensitivity in case of low loads. The C L uncertainties were estimated to be up to 8% for the SET1 and up to 3% for the SET2. For the drag, the uncertainty is up to 5% for both the wingsail configurations.

The comparison on the aerodynamic coefficients for the SET1 configuration ( Table 2 ) shows a discrepancy of 33.5% on C L between the experimental value and the The difference is even more elevated with a discrepancy of 44.9% on the drag coefficient. The wind tunnel numerical reproduction enhanced the CFD pre- dictions on the wingsail: the errors drop down to 9.3% on the lift and 26.9% on the drag. This enhancement in the numerical modeling can be further observed comparing the C p distribution over the three wingsail sections ( Fig. 9 ).

In the FSNS the suction peak on the upper surface of the airfoil is overestimated from a minimum of 0.3 to a maximum of 0.5 with respect to the experimental data (A 1 , B 1 and C 1 ). Furthermore the transition takes place quicker than in the experimental case (A 2 , B 2 , and C 2 ). These problems are completely solved in the WTNS.

Here the match with the wind tunnel data is good. The C p distribution on the pressure side is correctly reproduced as well as in the suction side. The suction peak has a maximum discrepancy of 5% (A 1 , B 1 , C 1 ) while the transition zone is delayed by 1 to 2% of the chord. On the trailing edge zone (A 3 , B 3 , C 3 ), both the numerical solutions keep a lower pressure value than the experimental case, where the pressure tends to increase in the last 10% of the local chord.

The FSNS cannot correctly reproduce the experimental conditions because of the confinement effects introduced by the wind tunnel walls. In Table 3 the flow deviation angle on the local xy plane is reported for both the numerical cases WTNS and FSNS and compared to the three reference sections, at y = 0 and at the upstream distance (-x direction) of 20% of the local chord from the local main L.E. The flow deviation for the freestream simulation is 5 °-7 °more elevated than in the wind tunnel case. The effect of the wind tunnel domain reduces the actual angle of attack felt by the airfoils. The suction pressure capabilities are then worsened by this incidence reduction affecting the lift capabilities of the entire wing.

In the SET2 case, the discrepancy on the drag coefficient is reduced from 30.4% to 16.4% using the wind tunnel modeling. The lift is underestimated by 10.4% in WTNS while is overestimated by 8.8% in FSNS ( Table 4 ). The C p distribution analysis on the three wing sections ( Fig. 10 ) is more complex. On the lowest section a correct match is found between WTT data and WTNS data, while the freestream case overestimates the C p (A 1 ). On this same section the transition point is correctly estimated with an error of only a few percent on the chord compared to the experimental case (A 2 ). For z * = 0.5, the pressure side distribution is correctly reproduced by the WTNS; on the suction side on the contrary the FSNS properly estimates the C p on the turbulent zone (B 2 ) while the peak suction zone (B 1 ) is overestimated by the WTNS and overestimated by the FSNS. On the highest section ( z * = 0.75), the best match with the experimental data on the suction side is obtained with the FSNS, while the WTNS underestimates the suction on the entire chord (C 1 , C 2 , C 3 ). The wind tunnel here also reduces the angle of attack felt by the wing ( Table 5 ).

The reason why WTNS does not perform better than FSNS across the entire wingspan is more closely investigated in the next section. The main reason is that the SET2 configuration has a particular flow pattern, with the flow attached only on the low half sections of the flap while on the mid-high sections the flow is separated. This different condition on the flap influences directly the jet flowing inside the slot dividing the main from the flap. In the attached case ( i.e. on the lowest section) the slot jet lies on the flap surface keeping the direction given to it by the geometry. In the separated case ( i.e. on the high sections), the jet is detached from the flap surface assuming a direction dependent on the size of the flap recirculation zone.

Because of the wind tunnel influence, the flow in the WTNS is less deviated in the slot, giving a jet that has a low tangential momentum component. On the contrary the FSNS, that has no straightening effects, predicts a jet with a higher tangential momentum component. The difference in jet deflection between the two cases was estimated to be 2 °.

As described by [START_REF] Smith | High-lift aerodynamics[END_REF] the jet direction modifies the T.E. condition on the main element, changing its circulation. A jet deflection enhances the main circulation and hence its lift. This explains the differences found in the two numerical cases but nonetheless it does not justify why the experimental data match better with the free-stream case. Another explanation comes from the wing deformation that occurs during the wind tunnel tests. Particularly, the upper part of the flap moves away from the main element, widening the slot size. Having a larger slot, the jet was characterized by a lower velocity, thus preventing its capabilities to deviate in the same way as in the WTNS, where the slot is narrower. In actuality, the flow jet is less deviated in the FSNS and in the WTT, improving the main circulation and therefore enhancing the pressure suction on the high wing sections. Despite this discrepancy, the WTNS is the most appropriate approach to reproduce the experimental case and the flow physics. Only the WTNS database remains to be further analyzed in the following sections.

Flowfield comparison between the numerical analysis in wind tunnel and the experimental data

The skin friction features of the wing flowfield have been compared between the WTNS and the viscous oil visualizations. The wingsail map comparison shows a qualitative view of the flow pattern on the suction side, especially regarding laminar to turbulent transition zones and boundary layer separations. The flow is laminar on the first part of the wing chord (the zone near the L.E.). The flow separates, creating a laminar bubble that extends until the flow reattaches on the wing surface after the transition in turbulent regime that has taken place. The oil visualization highlights the location of laminar bubbles and separation lines where the oil stagnates. In the numerical solution the laminar bubble and the separated zones have been detected by means of the skin friction coefficient. The laminar and turbulent zones were then detected with the intermittency factor ( Suluksna et al., 2009;[START_REF] Malan | Calibrating the γ -Re ϑ transition model for commercial CFD[END_REF].

For the SET1 configuration the flow map is reported in Fig. 11 , showing a good agreement on the main element between experimental and numerical fields. The laminar bubble position is well detected on the entire wingspan with a discrepancy by 2% to 11% of the chord c 1 ( Table 6 ). An exception exists at the wing root, where the 3D flow phenomena, due to the flow interaction between the wing and the disk interface, make it harder to predict the transitional region.

Table 6

Numerical/experimental comparison for the position and the length of laminar separation bubble on the main element for SET1 configuration. The laminar bubble length is overestimated by 5%-13% by the simulation. Previous studies already reported this behavior with the transitional model used [START_REF] Malan | Calibrating the γ -Re ϑ transition model for commercial CFD[END_REF][START_REF] Chapin | Aerodynamic study of a two-elements wingsail for high performance multihull yachts[END_REF]. After the transition in turbulent regime, the flow is attached all over the main surface except on the trailing edge region.

The transition model predicts a transition on the flap that is induced by a laminar separation at half of the chord. However, in the oil flow visualizations, the transition is detected near the flap L.E. with a thin laminar bubble. The transition model considers the value of the local turbulent kinetic energy, as predicted by the turbulence model, to estimate the value of the Reynolds number based on the momentum thickness. Indeed, the incorrect transition detection is related to the difficulty of the RANS-based numerical simulation to accurately estimate the turbulent kinetic energy at the frontier of the flap boundary layer, which is protected by the jet from the slot.

For the SET2 configuration ( Fig. 12 ), the numerical and the experimental data are in good agreement. The flow on the main has the same characteristics as in the SET1. The laminar bubble position is well detected on the entire wingspan with a difference of less than 4% of chord as reported in Table 7 . Its length is overestimated with a difference from 2% to 12%. The flow separates from the main at 90% of the chord.

The transition model well detects the flap laminar zone that lies between the 5% and 10% of the flap chord. The numerical simulation captures the different flow features along the wingspan. From root to half span, the flow is attached on most parts of the surface, with a separation that occurs at 95% of the flap chord. On the upper part of the flap, the flow is completely separated.

Table 7

Numerical/experimental comparison for the position and the length of laminar separation bubble on the main element for SET2 configuration. This behavior of the flow around the flap is linked to the size of the slot along the wingspan. The jet of the slot is sensitive to the gap size, modifying the location of the flow separation on the flap. These findings emphasized the slot effect in wingsail performances as it will be described in more details in the next section.

Jet slot analysis

In existing literature, the characterization of the slot is expressed by two parameters: the overlap ( o ) and the gap ( g ) distance. The first one expresses the horizontal distance between the main T.E. and the flap L.E.; it assumes negative values when the flap is placed rearward to the main T.E. The gap size is the minimum distance between the main surface and the flap surface. Instead of the gap size, the vertical distance between the main T.E and the flap L.E. ( y F ) can also be used as parameter ( Fig. 13 ).

The presence of the jet improves the high lift capabilities of the flapped wing configuration as reported by [START_REF] Smith | High-lift aerodynamics[END_REF] . The jet of the slot can be considered as a potential flow lying between the viscous shear layer of the main element and the flap boundary layer. The high velocity region of the jet, deviated by the flap geometry, helps to increase the circulation on the main element enhancing its lift. The jet influences the flow near the T.E. of the main, decreasing its pressure recovery demands, bringing down the possibility of a flow separation from the main surface. The interference with the T.E. of the main reduces the flow momentum nearby the flap L.E., dumping its pressure peak and hence the pressure gradient on the flap surface, delaying the separation of the boundary layer.

The flap is separated from the main by the slot, so its boundary layer is thinner compared to the case of a flap fully incorporated in the main wing. The stability of this boundary layer is improved by the "off-surface pressure recovery", more efficient than a conventional recovery in contact with the wall. This phenomenon is due to the interaction among the three flow layers on the upper flap surface. The entraining effect due to the viscosity causes a deceleration of the flow that takes place on a thick zone and not only in contact with the surface. The deceleration is hence less abrupt with consequent lower adverse pressure gradients.

Nevertheless these effects only exist if the slot is properly dimensioned. [START_REF] Smith | High-lift aerodynamics[END_REF] proposed a criterion for the slot dimension, i.e . the slot has to ensure that the wake of the upstream el- o/c 1 y F /c 1 20 °Woodward and Lean (1993) -1.25% to 0.25% 1.5%-3.25% 30 °Biber and Zumwalt (1993) -2% to -0.2% 1.75%-3.25% 40 °Woodward and Lean, (1993) -1.25% to 0.2% 0%-2.2% Biber (2005) -1.5% to -1% 0.25%-0.75% ement and the flap boundary layer do not merge. Otherwise the merging of these zones is responsible for a thick viscous layer that is more prone to separation. Once the separation has occurred on the flap surface, the performance of the high lift configuration is reduced, potentially leading to a massive stall. [START_REF] Biber | Hysteresis effects on a wind tunnel measurements of a two-element airfoil[END_REF] described a double stall behavior on GA(W)2 high lift configuration: the flow separates, at first, on the flap surface, causing a first loss in lift, but without influencing the flow on the main that will separate when the inlet flow angle is still increased. Before the first stall, the lift slope tends to increase [START_REF] Woodward | Where is the high lift today? A review of past UK research programs[END_REF]. In fact, the increase of the angle of attack leads to a thickening of the wake, pushing the jet of the slot on the flap surface and hence delaying the flow separation. Nevertheless the slope enhancement is characteristic of high flap deflection angles (30 °-40 °) and in a certain range of gap and overlap that defines the slot optimal size. A rearward movement of the flap from the optimum position, as well as an y F increase, leads to a large separation on the flap surface. A forward movement of the flap or a reduction in the vertical distance ( y F ) is not detrimental to the lift at low and intermediate incidences. Nevertheless at high incidences the lift slope increase does not take place, limiting thus the maximum lift coefficient.

The slope enhancement does not appear for the 20 °flap deflection angle [START_REF] Woodward | Where is the high lift today? A review of past UK research programs[END_REF]. The maximum lift gradually decays when the flap is moved away from this position. The optimal dimension of the slot gap depends on the Reynolds number. At low Reynolds number, the boundary layer thickens, so it reduces the slot size felt by the flow. So the optimum size must to be larger at low Reynolds number than at high Reynolds number [START_REF] Haines | Scale effects on aircraft and weapon aerodynamics[END_REF]. Furthermore on wingsails the slot dimension is dependent to the flap deflection angle, because of the flap kinematic. The slot dimension of the studied wingsail is represented through its y F and o distributions along the span for the two configurations, Fig. 14 . The normalized values ( y F /c 1 and o/c 1 ) are not constant on the wingspan because of the tipward chord reduction.

Clearly for the SET1 configuration, where the deflection angle is smaller, the distance from the main to the flap y F is small and thus well adapted to avoid the flow separation over the flap. However for the SET2 configuration, the slot is adapted only on the lowest section of the wing. On the upper part of the span, the distance from the main to the flap y F is important and corresponds to conditions where the boundary layer on the flap is separated [START_REF] Woodward | Where is the high lift today? A review of past UK research programs[END_REF]. As reported in Table 8 , for low flap deflection angles, the optimum overlap corresponds to negative or small positive values, while for high deflection angles, the optimum moves towards more negative values. The optimum distance y F is larger at low deflection angles than at high deflection angles. Because of its flap rotating mechanism, it is thus difficult to obtain the optimal value for all flap deflection angles. When the flap angle is increased, the overlap of the slot moves toward positive values and larger distance y F dimension, which is exactly opposite to what should be done to obtain the optimal size.

To investigate the flow in the slot gap region, the numerical simulation is compared with 2D PIV-based measurements. The investigated configuration is the SET2 configuration. The analysis of the flow is led at z * = 0.25 (attached boundary layer) and z * = 0.75 (separated boundary layer). The scalar map for the normalized xcomponent of the velocity U/U ∞ is presented for z * = 0.25 ( Fig. 15 ) and for z * = 0.75 ( Fig. 16 ). The velocity profiles, U = f(l) (with l the distance to the wall) and the turbulent kinetic energy, k = f(l) are plotted at 90% of the main chord and at 10% of the flap chord, at z * = 0.25 ( Fig. 17 ) and z * = 0.75 ( Fig. 18 ). The velocity is normalized with freestream velocity and the turbulent kinetic energy is normalized with the freestream kinetic energy. The distance from the wall surface ( l ) is normalized with the local boundary layer thickness ( δ BL ). The turbulent kinetic energy measured with 2D PIV only takes into account for axial and tangential fluctuating velocity components (so the spanwise component is not taken into account). To compare numerical predictions with measurements, turbulence kinetic energy must be scaled down by a factor of 2/3 (the turbulence model assumes turbulence is isotropic).

At z * = 0.25, close to the root where the boundary layer remains attached, experimental and numerical data are in good agreement, especially for the main element. On the flap, some discrepancies appear along the trajectory of the wake induced by the main element. Numerical simulation successfully predicts the jet but it fails to reproduce properly the mixing layer between the jet and the wake of the main element. The PIV scalar map shows a velocity deficit in the wake of the main element all along the flap chord. In the numerical simulation this deficit is observed only in the neighbors of the flap L.E. but it is then quickly dissipated. Since the wake merges with the flap boundary layer (in the numerical simulation), it makes the flow more sensitive to the adverse pressure gradients U/U ∞ -0.7 1.7 0 1 and hence to separation (as shown also in Fig. 11 ). The turbulent kinetic energy agrees well with the PIV and numerical solution as shown in Fig. 17. At z * = 0.75, close to the tip where the boundary layer is separated from the flap surface, experimental and numerical data are also in good agreement. The most important discrepancy is related to the prediction of the slot jet direction, which is oriented This deformation induces a negative overlap between the main element and the flap. It also reduces the distance y F between both elements. As previously discussed, this effect can delay the separation of the boundary layer on the flap. In Fig. 18 , experimental and numerical data are in good agreement, showing a separated boundary layer on the flap at 10% of the flap chord. However, the peak of turbulent kinetic energy, related to the mixing layer between the jet flow and the wake of the main element, is found closer to the flap wall in the case of numerical simulation. This observation confirms that the numerical simulation does not predict accurately the trajectory of the jet flow.

This discrepancy about the jet deviation does not depend on a fault of the numerical approach but rather on a modification of the geometry caused by the mock-up deformation during the experimental campaign.

Conclusions

The rigid wingsail is an effective propulsion system that enhances yacht performance. Such wingsails can operate in severe conditions, showing massive boundary layer separation, at moderate to high Reynolds number. The objective of this study was to better understand the flow physics of a rigid wingsail, at low and high flap deflection angles. A particular attention was paid to the behavior of the flow in the vicinity of the slot between the two elements. The investigations were supported by a wind tunnel campaign coupled with 3D unsteady RANS simulations on a scaled wingsail. Wind-tunnel measurements on the two-element wingsail, typical of an AC72 design, were performed at a Reynolds number Re = 3 × 10 5 (based on root chord). Measurements include: aerodynamic loads, steady pressure sensors, oil flow visualizations and PIV fields.

Two approaches were tested for the unsteady RANS simulations: the classical freestream conditions and the wind tunnel environment. Comparisons were carried out on these different situations emphasizing the following points:

1. The numerical predictions are improved when the wind tunnel environment is modeled (as generally reported for high lift configurations). Numerical simulations in freestream conditions overestimates aerodynamic coefficients (lift and drag) compared to experimental data. 2. The numerical simulations in wind tunnel environment demonstrate its capability to predict the attached and separated regions as well as laminar and turbulent regions. Comparisons with PIV measurements confirm the ability of unsteady RANS to predict the flow around wingsail both with low and high flap deflection angles. 3. Some differences have been identified on the prediction of the mixing layer between the jet flow and the wake of the main element. Such flow is known to be difficult for turbulence mod-els. Another source of discrepancy comes from the real geometry of the wingsail, which experience shape deformation during wind-tunnel tests. 4. The flow physics in the slot is a key element of wingsail performance.

Further numerical simulations are needed to enhance the prediction of the jet flow, with LES based methods that are able to deal with massively separated flows and mixing layers. The real slot geometry of the full-scaled wingsail, including its deformations, should be quantified to estimate the uncertainties associated to the real slot geometry. The effects of Reynolds number on the flow should be investigated, through numerical simulations on a full-scale wingsail ( Re = 0.53 × 10 6 in the present work while Re = [3 × 10 6 ; 10 × 10 6 ] for a full scale AC72).

Fig. 1 .Fig. 2 .

 12 Fig. 1. America's Cup catamaran AC72 propelled by wingsail foiling in the Bay of S. Francisco (photograph Carlo Borlenghi).

Fig. 3 .

 3 Fig. 3. Scheme of the S4 wind tunnel facility and main parameters.

Fig. 4 .

 4 Fig. 4. Box domain and mesh section for the freestream simulation.

Fig. 5 .

 5 Fig. 5. Lateral view of the wingsail and sections at half wingspan for Set1 (top) and Set 2 (bottom) configurations.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Section at y = 0 of the polyhedral mesh for the empty wind tunnel domain.

Fig. 8 .

 8 Fig. 8. Scalar map on the wingsail upper surface in the SET1 configuration, colored with the normalized distance to the wall y + .

Fig. 9 .

 9 Fig. 9. Comparison of the Cp distributions on the three sections of the wingsail for SET1.

Fig. 10 .

 10 Fig. 10. Comparison of the Cp distributions on the three sections of the wingsail in the SET2 configuration.

Fig. 11 .

 11 Fig. 11. Scheme of the different flow zones on the suction side of the wingsail in WTT (up) and in WTNS (down) for SET1.

Fig. 12 .

 12 Fig. 12. Scheme of the different flow zones on the suction side of the wingsail in WTT (up) and in WTNS (down) for SET2.

Fig. 13 .

 13 Fig. 13. Slot parameters: gap (g), overlap (o), y F .

Fig. 14 .

 14 Fig. 14. Distribution of the slot dimensions y F /c 1 and o/c 1 on the wingspan for the two configurations of the wingsail.

Fig. 15 .Fig. 16 .Fig. 17 .

 151617 Fig. 15. Scalar maps colored with the normalized velocity on SET 2 at z * = 25% from CFD (top) and PIV (bottom). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .

 18 Fig. 18. Comparison between CFD and PIV data at two stations on the main (up) and on the flap (down) for the SET2 configuration, at z * = 75%: velocity (left) and turbulent kinetic energy (right).

Fig. 19 .

 19 Fig. 19. Distribution of the slot dimensions y F /c 1 and o/c 1 on the wingspan in the ideal and in deformed case for the SET2 configuration.

Table 1

 1 Aerodynamic coefficients for SET2 configuration in WTNS at four different mesh refinements.

		MESH0	MESH1	MESH2	MESH3
	Cell count (Millions)	18	24	32	45
	C L	1 .239	1 .176	1 .124	1 .090
	C D	0 .150	0 .146	0 .143	0 .142

Table 2

 2 Aerodynamic coefficients comparison between experimental data and numerical solutions in both freestream and WT domains for SET 1 configuration.

		WTT	FSNS	WTNS	FSNS/WTT (%)	W TNS/W TT (%)
	C L	0 .773	1 .032	0 .845	+ 33 .5	+ 9 .3
	C D	0 .089	0 .049	0 .065	-44 .9	-26 .9

Table 3

 3 Flow deflection angles in the local xy plane at y = 0 and at the upstream distance of 20% of the local chord from the local main L.E. on the x-direction for the wing in WTNS and FSNS for the SET1 configuration.

					WTNS	FSNS
			z * = 0.25	7 °13 °z
			* = 0.50	15 °22 °z
			* = 0.75	13 °18 °Table
	4
	Aerodynamic coefficients comparison between experimental data and
	numerical solutions in both freestream and WT domains for SET 2
	configuration.
		WTT	FSNS	WTNS	FSNS/WTT (%)	W TNS/W TT (%)
	C L	1 .254	1 .365	1 .124	+ 8 .8		-10 .4
	C D	0 .171	0 .119	0 .143	-30 .4		-16 .4

Table 8

 8 Optimal values for y F and overlap sizes at different flap deflection angles.

Acknowledgments

This research was supported in part by ASSYSTEM France, which kindly provided PhD funding. Particular thanks to Laurent Neau, who made possible the creation of the PhD project on the wingsails.

Many thanks to the technical team of ISAE/DAEP, Patrick Chèze and Daniel Gagneux for the scale model manufacturing, Hervé Belloc, Valérie Allan, Henri Dedieu and Emmanuel Rivet for their valu-able support during the wind tunnel tests. Particular thanks to Philippe Barricau for his help during the PIV campaign.

The support provided by the computing center of the University of Toulouse (CALMIP) is also acknowledged (project p1425).