text:
PhD Camille Castells

François Richez

Michel Costes

Franc ¸ois Richez

Numerical Analysis of RPM effect on Dynamic Stall Phenomena on Helicopter Rotor at High Thrust Forward Flight

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Dynamic stall is one of the most limiting phenomena occurring during helicopter operation, especially in forward flight or in maneuvers. In these flight conditions, a dissymetry of airspeed is observed on the rotor disk due to the combination of rotating and advancing speed with respect to the azimuthal position of the blade. To balance the rotor loads and remove the corresponding rolling motion, the angle of attack has to be periodically modified. It is reduced on the advancing side where the airspeed is higher and it is increased on the retreating side where the velocity is lower. For the most severe flight conditions, the local angle of attack on blade sections exceeds the static stall angle and reaches a critical angle of attack. Beyond this angle, the boundary layer massively separates, the lift is considerably decreased and a strong negative pitching moment is observed. This phenomenon occurring at each rotor revolution causes vibrations, loss of the aircraft controls and eventually damages to the blade structure. It is called dynamic stall (Ref. 1).

Dynamic stall is still one of the most studied phenomena experienced by helicopter rotors since new rotorcraft are expected to operate with higher speed and thrust. However, it remains a challenging phenomenon of the rotorcraft environment for numerical prediction since it involves several mechanisms which are not well understood and challenging to predict. One can 1 name three-dimensional unsteady separated flow, transonic flow including shocks, reverse flow, blade-vortex interaction, vortex shedding or fluid-structure coupling.

A renewed interest on the problem could be noted over the last twenty years thanks to the progress of CFD in dealing with viscous separated flows, progressively switching from simplified cases to more realistic ones. Until very recently, most studies dedicated to dynamic stall were carried out on simplified configurations, typically twodimensional airfoils in oscillating pitching motion simulating the evolution of the angle of attack during a rotor revolution (Ref. 2). Experimental (Refs. 2, 3) and numerical (Refs. [START_REF] Pailhas | Experimental investigation of dynamic stall[END_REF][START_REF] Richter | Improved two-dimensional dynamic stall prediction with structured and hybrid numerical methods[END_REF] studies allowed scientists to investigate and obtain a better understanding of the different mechanisms occurring during a dynamic stall event. It highlighted the hysteresis aspect of the lift and the pitching moment. It also brought out that the sharp decrease in pitching moment is observed before the decrease in lift. It is explained by the shedding of the dynamic stall vortex travelling downstream from leading edge to trailing edge. Two-dimensional simulations allowed the determination of the numerical requirements needed to capture dynamic stall events for a reduced time-cost. One can mention the work of Refs. 5 and 6 on the time and space resolution, the work of Srinivasan et al. (Ref. [START_REF] Srinivasan | Evaluation of turbulence models for unsteady flows of an oscillating airfoil[END_REF]) on turbulence modelling or the investigations carried out on the influence of the prediction of boundary layer transition on stall (Refs. [START_REF] Costes | Numerical investigation of laminar/turbulent transition effects on the dynamic stall of an oscillating airfoil[END_REF][START_REF] Richter | Experimental investigation of unsteady transition on a pitching rotor blade airfoil[END_REF]. Reference 10 also highlighted that two-dimensional numerical simulations cause an earlier onset of the stall event, overpredict the separation and delay the reattachment with respect to similar three-dimensional configurations.

Then, three-dimensional configurations were investigated to take into account the influence of the finite span character of the blade and the impact of the tip vortex on stall onset. Le Pape et al. showed experimentally on a OA209 wing that three-dimensional mechanisms are present during a stall event (Ref. 11). The tip vortex of the wing produces a downward flow which tends to delay the stall onset. Thus, the flow separation occurs at the inner part of the wing and then spreads to the tip where the tip vortex delays flow separation. Numerical studies (Refs. [START_REF] Costes | Numerical investigation of three-dimensional effects during dynamic stall[END_REF][START_REF] Kaufmann | Numerical investigation of three-dimensional static and dynamic stall on a finite wing[END_REF][START_REF] Jain | High-resolution CFD predictions for static and dynamic stall of a finite-span OA209 wing[END_REF] In realistic rotor environment, the dynamic stall mechanisms are expected to be quite different from the isolated wing since blade-wake interaction or blade articulation and deformation can influence the stall onset by modifying locally the angle of attack. Bousman (Ref. 17 [START_REF] Ortun | Rotor loads prediction on the ONERA 7A Rotor using loose fluid/structure coupling[END_REF][START_REF] Richez | Numerical analysis of dynamic stall for different helicopter rotor flight conditions[END_REF][START_REF] Grubb | High fidelity CFD analyses of dynamic stall on a four-bladed fully articulated rotor system[END_REF][START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF] on different stalled flight conditions highlights the capability of numerical tools to predict dynamic stall on rotor configurations. They provide some clues on the several mechanisms involved in the stall onset such as blade-vortex interaction in a high angle of attack area (Refs. [START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF][START_REF] Richez | Analysis of dynamic stall mechanisms in helicopter rotor environment[END_REF] or the torsional response of the blade (Ref. [START_REF] Richez | Analysis of dynamic stall mechanisms in helicopter rotor environment[END_REF].

The objective of the present work is to get a better understanding of the mechanisms occurring during stall events in rotor environment, especially the influence of the RPM on stall onset and strength. Three flight conditions with potentially different stall mechanisms are compared using coupled CFD/CA simulations.

EXPERIMENTAL DATABASE

A wide range of flight conditions of the model-scale 7A rotor, represented in Fig. 1, have been tested in the transonic closed circuit S1MA wind tunnel at ONERA Modane test center, France, in 1991 (Ref. [START_REF] Crozier | Recent improvements in rotor testing capabilities in the ONERA S1MA wind tunnel[END_REF]). The 7A rotor, of radius R=2.1m, is composed of four articulated blades of constant chord c = 0.14 m and constant twist rate of -8.3 deg/R. The blade is made of two airfoils. The 13% relative thickness OA213 airfoil is used from inboard to 0.75R and the 9% relative thickness OA209, from 90% of the rotor radius to the tip. An extensive database is available including high thrust and high speed cases. It provides unsteady pressure measurements, integrated section normal forces and pitching moments at the five blade radial stations r/R = 0.5, 0.7, 0.825, 0.915, and 0.975 as well as blade bending moments and motion (rigid angles and elastic deformations).

Three configurations have been selected in the database for the numerical investigation presented in this work. They reflect the influence of the rotor angular speed on the dynamic stall mechanisms. All selected conditions have the same propulsive and lift forces. Table 1 summarizes the flow parameters for each flight condition investigated.

The decrease in RPM increases the dissymetry of velocity on the rotor disk. Thus, to maintain a trimmed rotor state, an increase of the angle of attack is necessary when the RPM is decreased. An increase of the pitch angle and of the flapping amplitude (Modane law : β 1c = θ 1s ) is actually obtained (Table 1). It results in an amplification of the stall event when the RPM is decreased. Dynamic stall is characterized by a severe loss of lift and a sharp negative pitching moment on the retreating blade, for 180 deg ≤ ψ ≤ 360 deg. Figure 2 ). One can also notice that for the most critical case (LRPM), a second strong stall event occurs at 330 deg for the r/R = 0.915 section which is not seen in the other cases. Experimental data are not available for the High RPM case at r/R = 0.7 section due to pressure sensors failure. Thus, this section will not be presented in the paper.

The objective of the present work is to get an understanding of the mechanisms occurring during these stall events by using numerical simulations. The HOST comprehensive analysis code from Airbus Helicopters (Ref. [START_REF] Benoit | HOST, a general helicopter simulation tool for germany and france[END_REF]) is used to take into account the blade motion and deformation in the coupled simulations. It is based on the lifting line theory for aerodynamics and beam theory with modal decomposition to describe the blade structure. Seven modes and seven harmonics have been used to characterize the blade deformation. The aerodynamic loads provided by the CFD solver are used by the comprehensive code to compute the trim control angles and also the kinematics and deformations of the blades, which are, then, used by elsA to prescribe the motion and deformation of the blade grids. The coupling between codes is made at each rotor revolution. The simulations are stopped when the rotor control angles and the airloads reach an equilibrium state. For the present computations, such a convergence typically requires 10 to 15 revolutions.

COMPUTATIONAL METHODOLOGY

VALIDATION

The converged trim control angles obtained by the coupling process are compared to the experimental data available for each flight condition. The numerical calculations are in good agreement with the wind tunnel results as shown in Table 3 The trim angles are well estimated by the coupling method with a difference of less than 0.3 deg with experiment for all control angles. The test-stand used during the experiment has been taken into account in the computation (Fig. 4) since Ortun (Ref. 20) shows that it improves the lateral cyclic pitch angle estimation of the simulation by deflecting the flow in the rotor disk and so modifying the local angle of attack. The test-stand creates an upwash flow on the front of the rotor, ψ = 180 deg, and a downwash flow on the rear part, ψ = 0 deg, thus inducing a decrease of the local angle of attack at ψ = 0 deg and an increase at ψ = 180 deg which is directly balanced by the lateral pitch angle.

The sectional normal force and pitching moment coefficients (respectively M 2 C n and M 2 C m ), which are available in the experimental database, are compared to the computed results in Figure 5 to 10. In order to avoid a systematic error to get these quantities from elsA, the same integration procedure applied in the tests was used by interpolating the computed pressure at the location of the unsteady pressure transducers.

For the moderate RPM test case, the dynamic stall event is captured by the simulation. The decrease of M 2 C n and M 2 C m are well predicted by the coupled calculations (Figs. 9 and6). The mimimum of M 2 C n , obtained at ψ = 300 deg for the r/R = 0.915 section, is slightly overestimated but the stall event is well predicted in amplitude and in phase. The pitching moment coefficient is more sensitive and more difficult to capture, but the agreement between simulations and experiments is satisfactory since the moment stall occurs at the right azimuthal position. The minimum of M 2 C m is reached at ψ = 290 deg for the r/R = 0.915 section. One can notice that the minimum value reached is underestimated (more negative) by the CFD/CA computation. A second small variation of M 2 Cm, which is visible but less intense in the experimental data, is predicted by simulation for r/R = 0.915, ψ = 330 deg.

It is probably due to the slight overestimation of the pitch angle and of the shaft.

For the high RPM case, the stall observed is weaker than the one observed in the MRPM case (Fig. 3). Figure 5 shows that the sectional normal force coefficients are well captured by the calculation. However, one can notice, in Figure 8, that the negative pitching moment peak induced by stall is almost missed by the simulation. A slight decrease in the pitching moment is indeed observed at stall position ψ = 290 deg for r/R = 0.915 but it is strongly underestimated with respect to the experimental data. This kind of behaviour has already been observed in two-dimensional configurations. Light stalls where the pitching angle oscillations remain close to the stall angle are more difficult to capture by numerical simulations.

The two dynamic stall events observed for a lower RPM are captured with a good agreement with experimental results (Figs. 7 -10). Both of the events are however overestimated (more severe) in terms of pitching moment coefficient but are well captured in terms of azimuthal phase of the stall onset, the first stall event starting at ψ = 240 deg and the second one at 330 deg at the section r/R = 0.915 (Fig. 10). Therefore it can be concluded that the loose coupling method correctly captures dynamic stall events for these three different flight configurations.

Besides the airloads, the coupled procedure shows a good capability to predict the structural loads of the rotor blades. is illustrated in Figure 15. They reveal four flow separation regions for all cases. They are circled in black. They also exhibit an increase of the global separated flow area on the rotor disk with the decrease in RPM. The reverse flow region is a region close to the retreating blade root where the flow is coming from the trailing edge to the leading edge. From the definition of the separation criterion based on the friction lines, the reverse flow is detected as a boundary layer separation but it will not be considered in the following discussion.

The first flow separation occurring on the blade is located on the inner part of the blade (0.4 ≤ r/R ≤ 0.8) in the third quadrant of the rotor disk, 180 ≤ ψ ≤ 270. This region is referred to as "A" in Figures 14 and15. The boundary layer separates at the trailing edge at r/R = 0.6, ψ = 150 deg and moves towards the leading edge when the blade pursues its revolution. Figure 14 also highlights that the flow separation moves further upstream when the RPM is decreased. Thus, the flow separation point does not reach further than the quarter of chord in the high RPM case, x sep /c = 0.3 at ψ = 250 deg. In contrast, it reaches the leading edge, x sep /c ≈ 0, for the two lower RPM cases and at an increasingly early azimuthal position of the blade. For the r/R = 0.5 section, the separation point reached the leading edge at ψ = 240 deg for the moderate RPM case and at ψ = 225 deg for the lowest RPM case. The separation also spreads radially, contaminating more external sections of the blade. Figure 15 shows that, for the two lowest RPMs cases, this region leads to a massive flow separation (l sep /c ≥ 0.9) which is responsible of the stall events observed in Fig. 5 -10, at r/R = 0.5, while no stall is observed in the HRPM case in this section.

A second flow separation region occurs at the blade tip for 210 deg ≤ ψ ≤ 310 deg and is referred to as "B" in Figures 14 and15. Firstly for ψ ≤ 270 deg, the separation suddenly appears on the leading edge of the blade and immediately reattaches as illustrated in Figure 15, the separation flow length being less than 5% of the blade chord. However, as the blade reaches a critical azimuthal position, the bubble bursts for the two lowest RPM cases. Figure 15 exhibits the massive flow separation (l sep /c ≈ 1) leading to the stall events seen on these two cases at the tip region of the blade in Fig. 5 -10. It occurs at ψ = 290 deg on MRPM and at ψ = 270 deg on LRPM.

The bubble does not burst in the high RPM case since the flow separation remains smaller than 10% of the blade chord. The next flow separation area, referred to as C, is located at the end of the fourth quadrant (330 deg ≤ ψ ≤ 360 deg) also at the blade tip region (Fig. 1415). For all the cases, the boundary layer first separates and immediately reattaches leading to a small separation bubble. As the azimuth increases, the bubble expands and gives a massive flow separation only for the LRPM test case. This flow separation leads to the second stall event observed on the airloads in Fig. 7 and 10. For the moderate RPM case, the separated flow reaches 30% of the blade chord which explains the small pitching moment fluctuations seen in Fig. 9, at r/R > 0.915 and ψ = 330 deg. For the HRPM case, the separation bubble of region B (Fig. 14) does not reach a length sufficient to influence the loads since l sep /c < 5% (Fig. 15).

The Flow Separation region D (Fig. 14 In this section, the effect of the BVI on stall event of region A is investigated for each rotor configuration. An analytical tool is implemented to estimate the location of the interaction between a blade and the tip vortex of the preceding blades. To do so, an analytical prediction of the position of the tip vortices emitted is computed. A helicoidal shape of the vortex driven by the rotational velocity of the rotor and the translational speed of the aircraft is assumed. The lead-lag and flap motions of the blade are considered in the computation as well as the shaft angle. The radial position of the interaction with the blade can then be easily determined. For clarification purposes, a snapshop of the analytical vortices for a blade at ψ = 280 deg for the MRPM case is shown in Figure 16. A good agreement between the isocontour of Q criterion computed from the elsA/HOST simulation and the analytical vortices can be noted. The intersection point between the blade and the tip vortex paths at the retreating blade is cercled in black in Fig. 16.

The radial position of the intersection of the blade with the helicoidal vortices has been computed for each azimuth. These intersection lines are represented with black lines overlayed on the rotormaps of the chordwise flow separation point in Figure 17 and on the rotormaps of the chordwise flow separation length in Figure 18. For all cases, it is clear that the blade-vortex interaction and the separation flow development are linked, especially in region A.

The interaction occurring in region A has started at the tip of the blade at the azimuth ψ = 90 deg, when the blade is set at a low angle of attack. Thus, the vortex passage does not induce a sufficient rise of angle of attack to trigger boundary layer separation. When the blade goes on with its revolution, the pitch angle is increased and the vortex influence becomes more critical. A flow separation begins to develop from the trailing edge at r/R = 0.6 at the interaction point for ψ ≈ 180 deg and moves towards the leading edge as ψ and the angle of attack increase. The lower the RPM is, the earlier the stall occurs since the kinematic-induced angle of attack is higher when the RPM is decreased. Stall onset in region A is marked by a + for the two lower RPM cases in Fig. 17 and 18. For the HRPM test case, the sign + shows the region where the flow separation length reaches its maximum. Figure 19 illustrates the isocontour of Q criterion for all the cases for a blade at the position marked by the + in Fig. 17 and 18. The dynamic stall vortex characterized by a high level of turbulent kinetic energy can be distinguished from tip vortices that are characterized by low level of turbulent kinetic energy. At such azimuthal position, the stall events observed in Fig. 5-10, at r/R = 0.5, start growing and can be identified in Fig. 19 by the high level of turbulent kinetic energy k. The stall develops at the blade-vortex interaction point and then, as ψ is increased, spreads radially around the BVI as illustrated in Fig. 17 and18.

Analysis of Flow Separation in Region B

Figures 17 and 18 also show that the region B occurs in an area where the preceding blade tip vortex is close to the blade. A small separation bubble at the leading edge is observed before the BVI, for ψ < 270 deg for all cases in Fig. 17 and 18, but this separation length remains smaller than 10% of the blade chord. When the tip vortex comes closer to the blade, it seems to trigger the bubble bursting for the MRPM case. In this case, the sudden growth of the bubble size is occurring at the time of the interaction between the blade and the tip vortex (ψ = 280 deg). In Figures 17 and18 the point marked by a • represents the stall onset. It occurs at the azimuth of the blade-vortex interaction for the MRPM case. However, for the LRPM case, the stall onset is occurring before the bladevortex interaction. In Figure 20, the isocontour of the Q criterion for a blade at the • position is represented for all the cases, ie. at the stall onset for the two lower RPM cases and at the larger separation bubble for the HRPM case. At this position of the blade, the tip vortex is close to the blade tip where the stall vortex is growing, especially in the two lower RPM cases. The stall vortex on the blade surface is recognizable by the high level of turbulence energy. For the LRPM case, as observed in Fig. 17 and 18, the stall is triggered when the tip vortex is still far from the blade. Thus, the stall onset cannot in this case obviously be attributed to the blade-vortex interaction.

In this region, the mechanisms at the origin of stall appears to be more difficult to analyze than in the inner part of the blade (region A). To get a better understanding of the correlation between the motion of the blade and the stall events, further analysis of the phenomena occurring at the tip region of the blade is necessary. The azimuthal evolution of the pitch angle, the Mach number and the torsion angle at the r/R = 0.915 section are represented in Figures 22,[START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF] and 24 respectively for the HRPM, MRPM and LRPM cases. The sectional pitching moment coefficient and chordwise flow separation areas are also illustrated in these figures in order to compare the azimuthal phase of the phenomena. In order to facilitate the reading, the flow separation regions B, C and D observed at the blade tip in Fig. 14 and 15 are grayed.

In Figure 24, the pitch angle is found to be higherfor the LRPM case than the two other cases (Fig. [START_REF] Grubb | High fidelity CFD analyses of dynamic stall on a four-bladed fully articulated rotor system[END_REF][START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF]. It reaches 18 deg at the stall onset (ψ = 260). It is higher than the pitch angle obtained at the intersection point between the blade and the tip vortex leading to stall for the MRPM test case, θ < 17 deg at ψ = 270. Thus, it can be assumed that the variation of the angle of attack induced by the tip vortex could be enough to trigger the stall in this case, even if the vortex is still distant from the blade.

In Fig. 22, 23 and 24, the stall onset in region B is clearly seen in MRPM and LRPM cases. The separation bubble bursts at ψ = 270 deg for MRPM and at ψ = 260 deg for LRPM leading to a strong negative pitching moment. This excites the blade torsion mode, as a local minimum of the torsion angle is reached directly after the stall event. For the HRPM (Fig. 22), the separation bubble does not burst, probably due to the lower pitch angle compared to the two other cases. Thus, it does not lead to stall. Therefore, the amplitude of the blade torsion variation is smaller than for the lower RPM cases.

Analysis of Flow Separation in Region C

The third flow separation area, C, located at the end of the fourth quadrant of the rotor disk, in Fig. 14 and 15 also concerns the blade tip section. The flow separation at r/R = 0.915 is represented in Figures 22,23 and 24. At such azimuthal position, the pitch angle is still high. The blade undergoes a flapping down motion due to the Modane law which imposes β 1c = θ 1s . Furthermore, the torsional deformation angle of the blade is at a local maximum . Thus, one can assume that both kinematic and deformation of the blade contribute to increase the angle of attack in region C. The process assumed through the rotormaps analysis 17 and 18 is very similar to the one observed in region B. A small leading edge separation bubble is predicted for all cases at the end of the fourth quadrant which bursts at the blade-vortex interaction for the lowest RPM cases. In Figure 21, the isocontours of the Q criterion for the blade at ψ = 330 deg, marked by a in Fig. 17 and 18, are presented for all the cases. This azimuthal position has been chosen to evaluate the mechanisms involved in the load variations seen at the blade tip r/R > 0.915 in Fig. 8, 9 and 10. A new blade-vortex interaction takes place when the blade is at ψ = 330 deg. Indeed, in this azimuthal range, the blade encounters the tip vortex emitted by the previous blade half a revolution earlier interacting with the previous dynamic stall event (B) which was convected behind the rotor. The vorticity in the vortex impacted is higher in the LRPM case (Vorticity ≈ 350 -400s -2 ) than in HRPM case (Vorticity ≈ 250 -300s -1 ). It probably contributes to increase the vortex-induced angle of attack. Thus, itcan be concluded that the mechanisms triggering flow separation in region C is a combination of several phenomena. The bladevortex interaction coupled with the positive torsion angle and the high pitch angle (Fig. 22, 23 and 24) probably leads to this second stall event, observed in the sectional airloads (Fig. 5 -10, r/R = 0.915 and r/R = 0.975).

Analysis of Flow Separation in Region D

Flow Separation D occurs in a high Mach number region. The local Mach number is increasing due to the helicopter speed which is added to the rotational velocity of the blade until it reaches its maximum at ψ = 90 deg. The evolution of the local Mach number and the torsion angle at r/R = 0.915 is presented in Fig. 22 for the HRPM case, in Fig. 23 for the MRPM case and in Fig. 24 for the LRPM case.

A flow separation is observed in this region for all cases. It occurs earlier with a high RPM because the shock appears earlier due to a high Mach number coupled with relatively high angles of attack. Figures 22,[START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF] and 24 also illustrate the evolution of the chordwise position of supersonic areas at r/R = 0.915 section for the three test cases. K p,crit is the pressure coefficient corresponding to a local Mach number equal to one. Thus, supersonic regions are characterized by K p -K p,crit < 0. Figures 22, 23 and 24 also include the location of the separation and the reattachment points of the different flow separations previously described. In region D, the small separation bubbles appear at the root of the shock wave, defined as the recompression edge of supersonic regions.

Two events are observed in region D for the LRPM case. Each area also coincides perfectly with a positive torsion peak seen in Figure 24. In the first event, the torsion adds 2 deg to the local angle of attack of the blade, significantly higher than for the two other RPM (Fig. [START_REF] Grubb | High fidelity CFD analyses of dynamic stall on a four-bladed fully articulated rotor system[END_REF][START_REF] Chaderjian | Navier-stokes simulation of UH-60A rotor wake interaction using adaptive mesh refinement[END_REF]. This high torsion angle is probably induced by the second stall event in region C which excites the blade torsion mode. It surely promotes the flow separation in region D.

Separation flows occuring in region D does not lead to stall in these three configurations since sectional normal forces and pitching moments do not show strong variation in these azimuthal and radial positions.

CONCLUSION

Loose coupling simulations between CFD code elsA and CA code HOST have been carried out for three helicopter flight conditions undergoing stall events. They correspond to the same flight condition of the helicopter at high speed and high thrust, with three different main rotor RPM regimes. Numerical parameters and turbulence model have been selected following the best practices determined in prior studies on pitching wings and helicopter rotors in stall conditions. For all cases, the numerical simulations have shown satisfactory agreement with the experimental data with regard to trim parameters (less than 0.3 deg difference with experiment for all control angles), airloads and structural loads (dynamic stall events correctly predicted in magnitude and phase). Different stall events have been distinguished for each case. The stall intensity increases with the decrease in the RPM. A trailing edge stall is observed for all cases in the inner part of the blade. For the two highest RPM cases, the blade tip undergoes one stall event at the beginning of the fourth quadrant of the rotor disk while it sees a double stall in the fourth quadrant for the lowest RPM case. The mechanisms that lead to these different stall events have been deeply investigated in this work. The onset of stall has been associated to the impact of the tip vortex generated by the other blades. The trailing edge stall and the strong leading edge stall occurring at the blade tip are triggered by the angle of attack induced by the vortex passing close to the blade. Stall only occurs when the kinematic induced angle of attack added to the vortex induced angle of attack is high enough. The second stall event in the low RPM case is connected to a second blade-vortex interaction. It takes place when the blade meets the tip vortex coming from the previous blade. Further investigations will be necessary to determine the influence of the vortex features on the dynamic stall onset. It also appears that the torsional deformation of the blade is sensitive to stall events occurring on the rotor. The dynamic stall, generating periodic strong negative pitching moments, induces negative torsion angles. Thus, it excites the torsion mode and leads to positive torsion angle peaks in the first quadrant. It was shown that, for these flight conditions, these local torsion maximums are correlated to flow separations and, thus, to load variations.

Author contact: Camille Castells camille.castells@onera.fr

  used Ref. 11 to validate the CFD methods and to provide the best set of numerical parameters and physical models for stall prediction. The next step in the dynamic stall studies is the inclusion of the blade rotation. Schwermer et al. (Ref. 15) established an experimental set-up to investigate a rotating blade in pitching oscillation motion. The effect of rotation is to keep the stall vortex closer to the blade surface with respect to a non rotating configuration which was also observed numerically by Gardner and Richter (Ref. 16).

  ) was the first to investigate the dynamic stall in complete rotor environment on three conditions extracted from flight tests of the UH-60A Black Hawk helicopter (Ref. 18): two maneuvers and one stalled forward flight. Potsdam et al. (Ref. 19) developed a fluid/structure coupling strategy between a Comprehensive Analysis (CA) code and a Computational Fluid Dynamics (CFD) solver. Three flight conditions from Ref. 18 including the stalled forward flight studied by Bousman were computed. The coupling computation gave a good agreement with respect to the experimental data available. Ortun et al. (Ref. 20) also validated the CFD/CA coupled simulation on 7A rotor configurations. Richez et al. (Ref. 21) provided the numerical requirements necessary to predict stall at high thrust forward flight. Recent works (Refs.

   depicts the evolution of the torsional moments obtained by the CA/CFD with respect to the experiment. Two sections (r/R = 0.3 and r/R = 0.65) of the blade are represented. The computations are in good agreement bwith the experiments. The 5/rev frequency content of the torsional moment is well predicted by the elsA/HOST computations for all cases. The amplitude are overestimated in the MRPM and the LRPM cases but the azimuthal phase correlation between the simulation and the experimental data is satisfactory for the two sections considered.ANALYSIS Description of the separated flow regionsIt was shown above that the numerical results are in good agreement with the experimental data in these three stalled configurations. Thus, it was concluded that the simulation correctly captures the mechanisms of the stall events and is representative of the phenomenon investigated. The numerical results are used to identify the stall regions and their features during these flight conditions. Rotormaps of the chordwise flow separation point position x sep and chordwise flow separation length l sep at each radial section and azimuthal step have been established. The separation point, x sep , is defined as the point where the tangential component of the friction vector becomes negative and, in the same way, the reattachment point is reached when it becomes positive again. The tangential unit vector is defined as shown in Fig.12. The leading edge is the origin of the local chordwise axis so x/c = 0 represents the leading edge (LE) position while x/c = 1 is the trailing edge (TE). Thus, x sep /c = 0 defines a flow separation at the very leading edge while x sep /c = 1 means that there is no flow separation. In the same way, the chordwise length l sep of the flow separation at each radial and azimuthal position is computed and l sep /c = 0 means that no flow separation is detected while l sep /c ≈ 1 defines a fully separated flow along the chord.Castillo et al. proposed, in Ref.[START_REF] Castillo | Separation criterion for turbulent boundary layers via similarity analysis[END_REF], a boundary layer separation criterion based on the shape factor: H sep = 2.76. This criterion has been validated on experimental results. However, this criterion is not able to detect clearly the reattachment and thus the separation flow length in these configurations. The criterion based on the friction lines, used to detect flow separation features has been validated by comparison with the one based on the shape factor used by Richez in Ref.[START_REF] Richez | Numerical analysis of dynamic stall for different helicopter rotor flight conditions[END_REF]. Figure13shows that the two different strategies of flow separation detections predict similar separation regions. The postprocessing tool used by Richez in Ref. 21 filtered some flow separation regions notably the small bubbles near the stagnation point, thus explaining the differences noticeable in Fig. 13. The rotormaps of the chordwise flow separation point position computed for the studied configurations are presented in Figure 14 and rotormaps of the chordwise flow separation length

Fig. 1 :

 1 Fig. 1: 7A rotor mounted in ONERA Modane S1MA Wind Tunnel.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Experimental sectional normal force coefficient M 2 C n .

Fig. 4 :

 4 Fig. 4: Structured mesh: Chimera backgroung grid, blade mesh, rotor and test stand surface mesh.

Fig. 5 :

 5 Fig. 5: Sectional normal force coefficient M 2 C n for the High RPM case.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Sectional normal force coefficient M 2 C n for the Moderate RPM case.

Fig. 8 :

 8 Fig. 8: Sectional pitching moment coefficient M 2 C m for the High RPM case.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Sectional pitching moment coefficient M 2 C m for the Moderate RPM case.

Fig. 11 :

 11 Fig. 11: Torsion Moment.

Fig. 12 : 12 ShapeFig. 13 :Fig. 14 :Fig. 15 :

 1212131415 Fig. 12: Schematic representation of the method used for separated flow detection.

Fig. 16 :

 16 Fig. 16: Analytical tip vortex and isocontour of Q-criterion colored by the turbulent kinetic energy (MRPM).

Fig. 17 :Fig. 18 :

 1718 Fig. 17: Rotormaps of the chordwise position of flow separation x sep /c. Black lines represent the analytical position of the blade-vortex interaction. +, • and correspond to the position of the maximum separation length respectively in region A, B and C.

Fig. 19 :

 19 Fig. 19: Isocontour of Q criterion colored by the turbulent kinetic energy for the blade at the position marked by + in Fig 17-18.

Fig. 20 :

 20 Fig. 20: Isocontour of Q criterion colored by the turbulent kinetic energy for the blade at the position marked by • in Fig 17-18.

Fig. 21 :

 21 Fig. 21: Isocontour of Q criterion colored by the turbulent kinetic energy for the blade at the position marked by in Fig 17-18.

Fig. 22 :

 22 Fig. 22: Representation of the evolution of the Mach number, torsion angle, the pitching moment coefficient, the flow separation region and the chordwise position of the shock at r/R = 0.915 section (HRPM).

Fig. 23 :

 23 Fig. 23: Representation of the evolution of the Mach number, torsion angle, the pitching moment coefficient, the flow separation region and the chordwise position of the shock at r/R = 0.915 section (MRPM).

Fig. 24 :

 24 Fig. 24: Representation of the evolution of the Mach number, torsion angle, the pitching moment coefficient, the flow separation region and the chordwise position of the shock at r/R = 0.915 section (LRPM).

  represents the phase-averaged evolution of the normal force

		Low RPM Moderate RPM High RPM
		(LRPM)	(MRPM)	(HRPM)
	Ω, rpm	964	1022	1063
	M tip	0.616	0.646	0.676
	µ	0.315	0.3	0.288
	C L /σ	0.09	0.1	0.11
	C X /σ	0.0044	0.0046	0.0049
	α t , deg	-6.75	-6.7	-6.45
	θ 0 , deg	9.7	8.4	7.2
	θ 1 c, deg	3.6	3.15	2.5
	θ 1 s, deg	-4.4	-3.5	-3.05
		Table 1: Selected Flight Conditions
	coefficient of the five instrumented sections. Lift stall events
	are clearly observed for the three cases at the outer sections,
	r/R > 0.8 for 270 deg ≤ ψ ≤ 330 deg. Figure 3 illustrates the
	pitching moment coefficient for the available sections. The
	strong negative pitching moment at stall (M 2 Cm steep nega-
	tive peak) can be seen on the retreating blade, around 270 deg
	at the r/R = 0.915 section, for all the selected flight condi-
	tions. The reduction of RPM leads to an earlier stall onset due
	to higher angles of attack. For the highest rotational speed
	(HRPM), stall appears at the r/R = 0.915 section only at an
	azimuthal position of 280 deg while it starts at 240 deg for the
	lowest RPM tested (LRPM

Table 2

 2 validated for the test case at the moderate RPM flight condition (MRPM). It gives a solution which is weakly dependent with respect to the time resolution during the stall event.

	This work uses a loose coupling strategy between Computa-
	tional Fluid Dynamics (CFD) and Comprehensive Analysis
	(CA), similarly to the work presented by Potsdam (Ref. 19).
	The CFD part of the coupling procedure is performed with
	the unsteady Reynolds-Averaged Navier-Stokes (URANS)
	finite volume elsA (ONERA-Airbus-Safran property) solver
	developed by ONERA (Ref. 26). A Chimera approach is
	used, with near body curvilinear structured grids around
	each blade and around the experimental test stand, and
	background Cartesian grids (Fig. 4). The blade near body
	grid is composed of 5.7 million points (

). The k-ω Kok SST turbulence model is chosen because it has shown a good capability to predict flow separation in previous work

(Ref. 27)

. The time derivative is discretized with a second order implicit Gear scheme. At each time step, the non-linear problem is solved by an iterative Newton process. The time step is set to ∆ψ = 0.225 deg with 30 sub-iterations in the Newton iterative process in order to get second-order time-accurracy (Ref.

[START_REF] Liggett | Temporal convergence criteria for time-accurate viscous simulations of separated flows[END_REF]

). The time step used in the following results has been

Table 2 :

 2 Mesh Specification

Table 3 :

 3 . Trim Control Angles Comparison for each flight case

		LRPM MRPM HRPM
	Experiment: α t , deg	-6.75	-6.7	-6.45
	elsA/HOST: α t , deg	-6.44	-6.4	-6.1
	Experiment: θ 0 , deg	9.7	8.4	7.2
	elsA/HOST: θ 0 , deg	9.7	8.6	7.2
	Experiment: θ 1 c, deg	3.6	3.15	2.5
	elsA/HOST: θ 1 c, deg	3.9	3.36	2.8
	Experiment: θ 1 s, deg	-4.4	-3.5	-3.05
	elsA/HOST: θ 1 s, deg	-4.4	-3.7	-3.01

  -15) appears on the advancing side (0 deg ≤ ψ ≤ 90 deg) in the tip region of the blade. This flow separation does not reach a size which can lead to stall. For all cases, Figure15shows that l sep /c < 10%c in region D. A second flow separation occurs in the region D in the LRPM test case. It appears at blade tip in the same region of the rotor disk, in the first quadrant.

	Analysis of Flow Separation in Region A
	In Reference 31, Droandi et al. showed experimentally and
	numerically that the incidence induced by the passage of a
	perpendicular vortex close to a NACA 23012 airfoil causes a
	flow separation onset leading to stall. Richez (Ref. 24) high-
	lighted that the blade-vortex interaction also triggers the dy-
	namic stall onset on realistic configurations of helicopter ro-
	tor. In the studied rotor configurations, the vortex is not al-
	ways perpendicular to the blade when the interaction takes
	place which is quite different from what was already investi-
	gated in two-dimensional configurations (Refs. 31, 32). The
	flight condition investigated in Ref. 24 corresponds the mod-
	erate RPM case considered in the present work.

A. Zanotti, G. Gibertini, and A. Mencarelli, "Experimental investigation of perpendicular vortex interaction over an oscillating airfoil in dynamic stall conditions," in 39th European Rotorcraft Forum, pp. 85-88, 2013.

ACKNOWLEDGMENTS

Part of this research was supported by a DGA scholarship. The studies presented in this article is making use of the elsA-ONERA software, whose the co-owners are Airbus, Safran, and ONERA


domain: spi phys

text:
S Dequand
email: sylvie.dequand@onera.fr

A Geeraert

G Mortchelewicz

INNOVATIVE TAIL CONFIGURATION: NUMERICAL VS WIND TUNNEL TEST DATA

Keywords: aeroelasticity, high-fidelity simulations, flutter, empennage

Aircrafts can be designed with different aft body configurations and all kinds of empennages exist. This paper focuses on the aeroelastic behavior of intersecting surfaces typically arising on aircraft tail, such as U-tail. It presents aerodynamic and flutter correlations between experimental and numerical results, in subsonic and transonic domains for a wind tunnel model representative of such aircraft tail. The aim is to evaluate the ability of our high-fidelity numerical tools to reproduce complex aerodynamic phenomena due to flow interactions between the different airfoil surfaces and to predict the sensitivity of the aeroelastic stability to the geometric parameters (dihedral, yaw angle). Steady and unsteady pressure coefficients along different slices of the airfoil surfaces have been measured and are compared to numerically predicted pressure coefficients. Results obtained at different aerodynamic conditions provide a good insight of the effect of flow interactions between the intersecting airfoil surfaces. The aeroelastic stability of the flutter model is also investigated and allows the assessment of different numerical approaches by comparing their ability to predict flutter mechanisms.

INTRODUCTION

An aircraft tail is composed of a vertical and a horizontal empennage with fixed and movable parts which enable to control the aircraft stability in yaw and pitch, and the flight symmetry. Empennages with two or three vertical stabilizers are for instance often used for multi-engine aircrafts; twin-boom aircrafts are generally used in fret; other tail configurations are the T-tail, U-tail, V-tail, cross-tail and low-tail configurations. Innovative tail configurations are more and more considered with the objective of reducing noise and fuel consumption and an accurate understanding of the complex aerodynamic phenomena involved in such geometries with intersecting surfaces is needed. The subject of flow interaction between aircraft parts has become of interest since the 50s with the investigation of steady aerodynamic of wingfuselage interaction [START_REF] Flax | The aerodynamics of low-aspect-ratio wings and wingbody combinations[END_REF], or VTP/HTP interaction ( [START_REF] Davies | Generalized aerodynamic forces on a T-tail oscillating in subsonic flow[END_REF], [START_REF] Stark | Aerodynamic forces on a combination of a wing and a fin oscillating in subsonic flow[END_REF]). The approaches for aeroelastic stability analysis are usually the Doublet Lattice Method (DLM) ( [START_REF] Albano | A Doublet Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows[END_REF], [START_REF] Murua | T-tail flutter: Potential-flow modelling, experimental validation and flight tests[END_REF]) and also CFDbased aeroelastic simulations ( [START_REF] Geuzaine | Aeroelastic dynamic analysis of a full F-16 configuration for various flight conditions[END_REF], [START_REF] Nikbay | Coupled analytical sensitivity analysis and optimization of three-dimensional nonlinear aeroelastic systems[END_REF]). In this paper, the study of innovative aft body configurations relies on a rich and thorough experimental database, allowing detailed data analysis of aeroelastic phenomena. Preliminary comparisons between measurements on the U-tail configuration (Figure 1) and numerical results have already been presented in [START_REF] Geeraert | Wind tunnel flutter tests of a U-tail configuration. Part 1: Model design and testing[END_REF] and [START_REF] Mamelle | Wind tunnel flutter tests of a U-tail configuration. Part 2: Experimental and numerical results[END_REF]. This paper will present aerodynamic and flutter correlations between experimental and numerical results for a wind tunnel model representative of an aircraft U-tail. The aim of this study is to assess different numerical methods for the flutter analysis of tail configurations, in subsonic and transonic domains. The first section of this paper will present a brief description of the Wind Tunnel Tests and the different geometries studied. Then, the different numerical methods used for the flutter analysis will be presented and finally, the results (pressure coefficients, flutter diagrams) will be shown in the last section.

WIND TUNNEL TESTS

The tail flutter model tests campaign with various geometric parameters and different aerodynamic conditions was performed in 2016 in the ONERA transonic pressurized wind tunnel S2MA, in the framework of the Cleansky research program, and thanks to a cooperation between RUAG, Dassault-Aviation and Onera. These Wind Tunnel Tests have received funding from the European Union's Seventh Framework Program (FP7/2007-2013) for the Clean Sky Joint Technology Initiative under grant agreement CSJU-GAM-SFWA-2008-001. The flutter model was a bizjet type configuration, equipped with a large number of pressure sensors and accelerometers and was specially designed to investigate flutter mechanisms in subsonic and transonic domains. It was a 1/8 scale model wall-mounted in the wind-tunnel test section and consisted of a fixed fuselage and a half empennage geometry, as the U-tail shown in Figure 1. Different empennage geometries were investigated by varying dihedral and yaw angles of the vertical part. Furthermore, different kinds of measurements were performed during the test campaign. A pitch motion could be applied to the model by means of a remotely controlled electrohydraulic system. This so-called "pressure configuration" yielded to steady and unsteady aerodynamic data in order to calibrate numerical data. For a given Mach number, the critical pressure was determined by continuously increasing the stagnation pressure. Several stabilized runs with stepping pressure were also recorded and analyzed to build the flutter curves. The Frequency Response Functions of the accelerometer responses during the stabilized runs were calculated relatively to the excitation actuator sensor and modal parameters were identified using the LMS Polymax software. The flutter model was designed in such a way that the flutter mechanism was the result of a coupling between two structural modes, namely the first bending and torsion modes of the HTP.

NUMERICAL APPROACH

Low-Fidelity coupling method

The flutter equation ( 1) is usually written in the Laplace domain using the Laplace variable p:

[𝑝 2 M + 𝑝C + K -𝑞 𝑑𝑦𝑛 GAF(𝑀, 𝑘)]X(𝑝) = 0 (1)
where M, K, C are respectively the modal mass, stiffness and damping matrices, GAF is the matrix of Generalized Aerodynamic Forces, 𝑞 𝑑𝑦𝑛 is the freestream dynamic pressure, X is the modal coordinate vector, M is the Mach number and k is the reduced frequency. The usual approach in aircraft flutter prediction is to compute the Generalized Aerodynamic Forces (GAF) in the frequency domain and then to solve the aeroelastic stability problem using the classical double scanning "p-k" method ([10], [START_REF] Liauzun | Method of Fluid-Structure Coupling in Frequency and Time Domains using Linearized Aerodynamics for Turbomachinery[END_REF]). Most of the time, linear methods in frequency domain, such as the Doublet Lattice Method (DLM) [START_REF] Albano | A Doublet Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows[END_REF], allow to compute the GAF and are valid in subsonic domain. This method is used in this paper and p rovides, additionally to the experimental database, an efficient comparison tool to validate our high-fidelity numerical methods, at least for low Mach numbers. The aerodynamic mesh used in the DLM and shown in Figure 2, is made of different panels representing the different surfaces of the flutter model. The grey part corresponds to the fuselage (immobile part) and is made of a horizontal plane with 46 x 5 panels and a vertical plane of 46 x 10 panels. The blue part represents the HTP and is made of 21 chordwise panels and 13 spanwise panels. The green part corresponds to the VTP and is made of 15 chordwise panels and 8 spanwise panels. The structural mode shapes are projected onto these two parts by means of a finite volume fitting method. The mode shapes on the grey part are forced to zero by a locally smoothing technique applied at the junction between the grey and blue parts.

High-fidelity method

High-fidelity simulations are performed with the elsA CFD code (Onera-Airbus-Safran property) [START_REF] Cambier | The Onera elsA CFD software: input from research and feedback from industry[END_REF], using its aeroelastic module "elsA-Ael". The aeroelastic computations are based on a modal approach for the structure. The URANS computations are performed with a centered finite-volume discretization on a structured mesh, a backward-Euler scheme with implicit LUSSOR relaxation for time integration, and a 2 nd order central Jameson scheme for spatial discretization. An aerodynamic mesh was built for each studied configuration. The meshes used for the highfidelity fluid-structure coupling simulations are structured meshes of about 8 millions mesh cells which are cut into a large number of blocks in order to use a maximum of processors and to reduce the CPU computation time. Figure 7 shows the wall surfaces of the U-tail configuration. All the URANS computations have been performed with the Spalart-Allmaras turbulence model, for Mach numbers of 0.7 and 0.8. The walls of the wind tunnel were not modelized and only the wall on which the model was mounted was taken into account by applying a wall boundary condition on the grey boundary shown in Figure 7. All the other farfield mesh boundaries were set to non-reflective boundary conditions.

Figure 3: Aerodynamic structured mesh used in the Hi-Fi numerical simulations

Forced Harmonic Excitation method (FHEM)

GAF can be computed by evaluating the unsteady aerodynamic response to any admissible motion of the structure, considered to behave linearly, in the frequency range of interest of the problem. In this Forced Harmonic Excitation method, several numerical simulations are performed, for a set of structural modes and for various values of excitation frequencies. While very robust, this harmonic motion approach requires however one unsteady simulation for each structural mode and for each frequency selected in the frequency range of interest, and is therefore very time consuming.

Pulse method

In the so-called "Pulse" method ([13], [START_REF] Silva | Identification of linear and nonlinear aerodynamic impulse responses using digital filter techniques[END_REF], [START_REF] Silva | Development of unsteady aerodynamic state space models from CFD-based pulse responses[END_REF]), the forced motion is not harmonic in time anymore but corresponds to a Heaviside function (step function from 0 to 1, at a given initial instant, Figure 4), applied to a structural mode. With such an approach, the impulse response of the system is computed in a single simulation for a full frequency range. The GAF can be written as a function of a transitory part 𝑇(𝑖, 𝑡) and the impulse response ℎ(𝑖, 𝑗, 𝑡):

𝐺𝐴𝐹(𝑖, 𝑗, 𝑛∆𝑡) = 𝑇(𝑖, 𝑛∆𝑡) + 𝑒 0 ∆𝑡 ∑ ℎ(𝑖, 𝑗, 𝑘∆𝑡) 𝑛 𝑘=0 ( 2
)
where ∆𝑡 is the time step, 𝑛∆𝑡 is the current time, and 𝑒 0 is an amplification factor.

The GAF can then be determined by means of a classical Fourier analysis of the time responses:

𝑓 ̂(𝜔) ≈ ∆𝑡 ∑ 𝑓(𝑘∆𝑡)𝑒 -𝑖𝜔𝑘∆𝑡 𝑛 𝑘=0 (3)
In this paper, three different methods have been therefore used to compute the GAF of the flutter model: high-fidelity aeroelastic simulations based on the FHEM and the Pulse method, and the DLM. Figure 5 shows a comparison of the modulus and the phase of the GAF computed with these three methods for a given flutter model configuration. The GAF modulus obtained with the Pulse method and the FHEM are very similar but differ from the DLM results. A very good correlation is however obtained for the GAF phase and shows that the three methods are predicting the same aerodynamic behavior.

Direct Coupling method (DCM)

An alternative approach to the previous presented methods (FHEM and Pulse method) is the so-called "Direct Coupling" method (DCM) in which several dynamic fluid-structure coupling simulations are performed at different stagnation pressure Pi, with an initial modal velocity applied to each structural mode. For the dynamic fluid-structure coupling simulations, the numerical technique consists in using an iterative process which is converged at each physical time-step. As schematically represented in Figure 6, the different steps performed at each time step are: 1. Computation of the flow solution (using a dual time-stepping technique) 2. Computation of the generalized aerodynamic loads (GAF) 3. Computation of the generalized coordinates, solving the mechanical modal system and using the implicit Newmark scheme for time discretization 4. Computation of the displacements of the aerodynamic surface grid, from the modal shapes and the generalized coordinates 5. Deformation of the aerodynamic grid The generalized aerodynamic forces (GAF) are deduced from the aerodynamic flow field and computed by integration of the pressure field, directly on the aerodynamic wall grid since the structural mode shapes are projected to the aerodynamic surfaces in a pre-processing step by using interpolation or fitting techniques ( ,,

S i A i  ΦΦ ): dS n p gaf Wall i , A i       (4) where i , A   , i=1
,N are the N first mode shapes, mapped on the aerodynamic surface grid.

The generalized coordinates q i are then computed from the dynamic modal system: 2 2

( 1, )

i i i i i i i i i i m q m q m q gaf i N        (5)
where m i is the generalized mass, ω i the pulsation, and  i the structural damping coefficient.

The wall displacements and velocities, A X  , and A X  are provided directly on the aerodynamic surface grid by: , 1 ( , ) ( ) ( )

N A i A i i X M t q t M      , , 1 ( , ) ( ) ( )
N A i A i i X M t q t M      (6)
For the deformation of the aerodynamic grid, an individual "modal mesh deformation" is first computed for each mode at the start of the coupled simulation, using a mixed analytical/transfinite interpolation technique. The instantaneous mesh deformation n A X  at time level n is then computed as a linear combination of these individual "modal" mesh deformations:

modal , 1 N n n A i A i i X q X      (7)
where n i q are the generalized coordinates solution of the dynamic modal system.

At last, a coupling iteration loop between the two disciplines has been setup within the physical time loop in order to get a converged coupled fluid-structure state at each physical time step. In practice, it is observed that two or three iterations are enough for most aeroelastic problems.

For each Pi, the time responses obtained with the DCM are assumed to be written as:

𝑠(𝑡) = 𝐶 + ∑ 𝑠 𝑘 (𝑡) = 𝐶 + ∑ 𝐴 𝑘 𝑠𝑖𝑛(𝜔 𝑘 𝑡 + 𝜙 𝑘 ) 𝑁 𝑘=1 𝑁 𝑘=1 𝑒 -𝛼 𝑘 𝜔 𝑘 𝑡 (8)
The set of parameters {𝐶, 𝐴 𝑘 , 𝜔 𝑘 , 𝛼 𝑘 , 𝜙 𝑘 } 𝑘=1 𝑁 is estimated by a fitting procedure based on the Levenberg-Marquardt algorithm [START_REF] Fletcher | Practical Methods of Optimization[END_REF]. Figure 7 gives an example of a typical time response obtained and its Fourier signal processing. This Direct Coupling method enables to determine the critical stagnation pressure at which the system becomes unstable, but requires however several computations and can become very time consuming.

RESULTS

Cp comparisons

Steady pressure coefficients

The steady pressure coefficients have been measured at different slices over the horizontal and vertical tail planes (HTP and VTP) (Figure 8). These experimental data can be directly used to calibrate our numerical data. Static fluid-structure coupling URANS simulations have been performed using our highfidelity CFD code elsA, in order to compare the numerical pressure coefficients to the measured steady pressure coefficients. Figure 9 shows the results obtained at Mach = 0.7 and Pi = 1.0 bar with and without taking into account the structural deformation. The curves are respectively denoted 'Static Coupling' and 'Rigid' in the different graphs representing different slices over the HTP and the VTP of the model. Contrary to the experimental steady pressure coefficients, the 'Static Coupling' pressure coefficients are different at the upper and lower HTP surfaces and show that the empennage numerically computed static deformation has a non-zero incidence. The comparison of the numerical static deformed empennage with the MDM data (Model Deformation Measurement) confirms indeed that the static deformation is not well predicted by our static aeroelastic simulations. This is due to the modal approach used, in which the structure is represented only by its first eleven modes. In a best practice, the static deformation should be computed by means of a flexibility matrix approach. However, this method requires a Finite Element Model of the structure, which was not available for this study. The MDM data show that the static deformation of the empennage is negligible (displacements lower than 0.5mm) and when considering the numerically predicted 'Rigid' pressure coefficients, a rather good correlation with the experimental data is indeed obtained at the different HTP and VTP slices. Figure 10 shows the effect of a higher Mach number on the steady pressure coefficients obtained at the 'H-A' HTP slice, which corresponds to the slice located closer to the fuselage. The main difference is that the pressure is not well predicted from the last 20% of the HTP chord. This is probably due to a problem of leakages which occurred during the wind tunnel tests, at the junction between the fuselage and the HTP. In order to take into account these leakages, a mass flow injection boundary condition is imposed at the first line of cells located on the fuselage around the HTP junction (thin white zone in the bottom right hand-side graph in Figure 11). This injection condition is defined as a blowing flow condition corresponding about to 5% of the flow velocity, towards the Ydirection. Figure 11 shows indeed that the injection boundary condition improves the prediction of the steady pressure coefficients. The next figures show the effect of a change in the empennage geometry on the pressure field.

The geometrical parameters considered are the yaw angle and the dihedral. Figure 12 shows that a positive yaw angle changes completely the pressure field distribution over the VTP surface, with the apparition of a recompression, probably due to the VTP incidence. The negative dihedral increases this phenomenon.

Figure 13 shows the steady pressure coefficients at different HTP and VTP slices for positive and negative yaw angles. When the VTP is directed towards the fuselage (positive yaw angle), a shock on the VTP occurs. As shown in Figure 14, the effect of the dihedral is less spectacular. The steady pressure coefficients are very similar, with only a slight increase when the dihedral is negative (towards the fuselage).

Unsteady pressure coefficients

The unsteady pressure coefficients have been measured by exciting the flutter model with a time-harmonic pitching motion. Different experimental pressure tests were made at frequencies ranging from 5Hz to 60 Hz, and amplitudes 0.1° and 0.2°, but only the 5Hz / 0.2° pitch motion is considered in the following results. Dynamic fluid-structure coupling simulations and particularly, harmonic forced motion simulations, have been performed by imposing a pitching movement to the HTP part. The Frequency Response Functions are computed relatively to the excitation signal and the real and imaginary parts of the pressure can be compared to the experimental data. The unsteady pressure coefficients obtained at different HTP and VTP slices are presented in Figure 15, for a Mach number of 0.7. A rather good correlation is obtained for the real part of the pressure coefficients. However, some discrepancies occur for the imaginary part of the pressure which is not well predicted around the HTP / VTP junction (V-E section). These discrepancies are even more important when increasing the Mach number, as shown in Figure 16, for a Mach of 0.8. The real part of the pressure coefficients is well predicted, except again near the HTP and VTP junction (V-E slice). An attempt of improving the results by adding an injection boundary condition (to take into account the leakages appeared during the wind tunnel tests), is not efficient.

Flutter diagrams

The different numerical methods previously presented are used to plot the evolution of the modal frequencies and damping ratios with the stagnation pressure. The resulting flutter diagrams are presented in Figure 17 and Figure 18 for the first two structural modes (first bending and first torsion modes of the HTP) for Mach numbers of 0.7 and 0.8. The results obtained correspond to the flutter model with a yaw angle of -2° and a dihedral of 0°, and the numerical results are compared to the corresponding experimental data.

In Figure 17, the DLM (blue lines) gives a very good correlation with the experimental data and predicts the flutter pressure with only an error of 4%. The Pulse method and the FHEM (green and pink curves) give very similar results with a flutter pressure underestimated by 9%, and show both a very good prediction of the modal frequencies. Then, even if the Direct Coupling method (DCM) (red lines) gives a very good prediction of the flutter pressure (with only an error of 2%), behavior of the modal frequencies and the damping is not correlated with the other numerical methods and with the experimental data. In Figure 18, the increase in Mach number yields to a drop in the flutter dynamic pressure.

For this Mach number of 0.8, all the numerical methods predict approximately the same value of the flutter pressure with an underestimation of about 15%. Once again, the Pulse method is able to predict very accurately the evolution of the modal frequencies, contrarily to the DCM. When taking into account a correction for the leakages appeared during the WTT measurements (with an injection boundary condition), the prediction of the flutter pressure can be improved only by 3%.

Mach effect

Figure 19 shows the evolution of the flutter pressure with the Mach number for the flutter model with a yaw angle of -2° and a dihedral of 0°. Once again, the DLM gives relatively good correlations with experimental data with an underestimated flutter pressure of from 5 to 15 %, even for high Mach numbers. The high-fidelity numerical tools (DCM and Pulse method) are not able to improve the DLM predictions.

Yaw angle effect

When the yaw angle is set to +2°, i.e. when the VTP leading edge is directed towards the fuselage, the flutter pressure is increased up to 2.0 bars. For this configuration, the comparison of our numerical data to the experimental data shows again some important discrepancies. The flutter diagrams in Figure 20 show that the DLM and the Pulse method are not able anymore to predict the good evolution of the modal frequencies and damping ratios (with an error of prediction higher than 30%). The DCM gives relatively better results (flutter pressure underestimated with an error of 15%) but still does not predict the good flutter behavior.

Dihedral effect

As for the positive yaw angle, a negative dihedral increases the flutter pressure because it produces an antisymmetric incidence on the VTP and then yields to a stabilizing moment. For this configuration with a dihedral of -10°, the flutter diagrams in Figure 21 show that the DLM and the Pulse method have some difficulties to predict the good flutter behavior, in particular they predict both a bad evolution of the modal frequencies. The damping ratios, and then the critical flutter pressure, seem however to be relatively well predicted.

CONCLUSIONS

The aim of this study was to compare and assess different numerical methods for the flutter analysis of tail configurations, in subsonic and transonic domains. Numerical results obtained with high-fidelity fluid-structure coupling methods such as the Forced Harmonic Excitation, the Pulse and the Direct Coupling methods are compared to Wind Tunnel Test data and to results obtained with the low-fidelity DLM method. Steady and unsteady pressure coefficients are first analyzed and enable to optionally correct our numerical data by adding an injection condition at the junction between the fuselage and the HTP, where some leakages had occurred during the WTT. This correction is particularly applied at a Mach number of 0.8 and slightly improves the prediction of the flutter pressure. Aerodynamic and flutter correlations between experimental and numerical results are quite good for a tail flutter model involving intersecting surfaces, with a yaw angle of -2° and a dihedral of 0°. However, for a positive yaw angle of +2°, the numerical methods are not able to correctly predict the flutter speed. In the same way, for a negative dihedral, the flutter behavior and particularly the evolution of the modal frequencies is not well predicted. The comparison between the different numerical methods yields to mixed results:

-the results obtained with the DCM are very disappointing for the flutter behavior prediction, whatever the aerodynamic conditions or the flutter model geometry. -the results obtained with the Pulse method are very promising and show that this method can efficiently replace the more classical, but very time consuming, FHEM. -the results obtained with the DLM are relatively very good even at high Mach numbers (Mach = 0.8), except for geometries with more complex flow behavior (positive yaw angle, negative dihedral), and show the robustness of this method. Further work is still required to complete the analysis of the WTT database, and in particular to investigate other flutter model geometries. This could help to better analyze and understand the reasons for the poor results provided by DCM. It would be further interesting to investigate the use of the Quadratic Constitutive Relation (QCR) instead of the Boussinesq relation in the Spalart-Almaras turbulence model in order to better predict the corner flow behavior near the intersecting surfaces.

Figure 1 :

 1 Figure 1: Flutter model wall-mounted in the WTT section. Scheme with notations and main geometric parameters considered in this paper.

Figure 2 :

 2 Figure 2: Aerodynamic panel mesh built for the DLM.

Figure 4 :

 4 Figure 4: Step function used as excitation signal in the "Pulse" method.

Figure 5 :

 5 Figure 5: Mach = 0.7, Pi = 1.0 bar. GAF evolution with the frequency computed with three different methods: DLM, FHEM and Pulse. (Yaw angle = -2° / 0° dihedral)

Figure 6 :

 6 Figure 6: Aeroelastic optional subsystem of elsA for aeroelastic simulations

Figure 7 :

 7 Figure 7: Examples of a time response and its Fourier signal processing.

Figure 8 :

 8 Figure 8: Location and notation of the different HTP and VTP slices

Figure 9 :

 9 Figure 9: Experimental and numerical steady pressure coefficients on different HTP and VTP slices at Mach=0.7 and Pi = 1.0 bar (Yaw angle = -2° / 0° dihedral)

Figure 10 :

 10 Figure 10: Mach effect on experimental and numerical steady pressure coefficients on the 'H-A' HTP slice (Yaw angle = -2° / 0° dihedral).

Figure 11 :

 11 Figure 11: Effect of an injection condition at the junction between the fuselage and the HTP (Yaw angle = -2° / 0° dihedral). Pressure coefficients on the 'H-A' HTP slice and pressure field near the injection zone.

Figure 12 :

 12 Figure 12: Effect of the yaw and dihedral angles on the steady pressure field at the cross section between HTP and VTP (Mach = 0.7 and Pi = 1.0 bar).

Figure 13 :

 13 Figure 13: Yaw angle effect on experimental and numerical steady pressure coefficients on different HTP and VTP slices at Mach=0.7 and Pi = 1.0 bar with a fixed dihedral of 0°.

Figure 14 :

 14 Figure 14: Dihedral effect on experimental and numerical steady pressure coefficients on different HTP and VTP slices at Mach=0.7 and Pi = 1.0 bar with a fixed yaw angle of +2°.

Figure 15 :

 15 Figure 15: MACH = 0.7 -Experimental and numerical unsteady pressure coefficients on different HTP and VTP slices obtained with a harmonic excitation signal of 5Hz and 0.2° (Yaw angle = -2° / 0° dihedral).

Figure 16 :

 16 Figure 16: MACH = 0.8 -Experimental and numerical unsteady pressure coefficients on different HTP and VTP slices obtained with a harmonic excitation signal of 5Hz and 0.2°. With and without the injection correction (Yaw angle = -2° / 0° dihedral)

Figure 17 :

 17 Figure 17: MACH = 0.7 -Flutter diagram: Comparison of the experimental and DLM data with different highfidelity numerical methods (Direct Coupling, Forced Harmonic Excitation, Pulse methods) (Yaw angle = -2° / 0° dihedral).

Figure 18 :

 18 Figure 18: MACH = 0.8 -Flutter diagram: Comparison of the experimental and DLM data with different highfidelity numerical methods (Direct Coupling, Forced Harmonic Excitation, Pulse methods). With and without injection correction. (Yaw angle = -2° / 0° dihedral).

Figure 19 :

 19 Figure 19: Evolution of the critical pressure with the Mach number. Comparison of the available numerical methods with experimental data. (Yaw angle = -2° / 0° dihedral).

Figure 20 :

 20 Figure 20: MACH = 0.7 -Flutter diagram : Comparison of the experimental and DLM data with different highfidelity numerical methods (Direct Coupling, Pulse methods). Without injection correction. (Yaw angle = +2° / 0° dihedral).

Figure 21 :

 21 Figure 21: MACH = 0.7 -Flutter diagram : Comparison of the experimental and DLM data with the Pulse method). Without injection correction (Yaw angle = +2° / -10° dihedral).

ACKNOWLEDGMENTS

This project has received funding from Cleansky-2 Joint Undertaking under the European Union's Horizon 2020 research and innovation programme under grant agreement N° CS2-AIR-GAM-2014-2015-01". This paper only reflects the author's view and the CS2 JU is not responsible for any use that may be made of the information it contains. More specifically, the authors would like to thank Gabriel Broux, Elsa Breus and Nicolas Forestier from Dassault Aviation for their help and fruitful discussions.