Thomas Bauvent
email: thomas.bauvent@cea.fr

Paola Trotti

Olivier Billoint

Jean-François Nodin

Yasser Moursy

Gabriel Molas

Gaël Pillonnet

J-F Nodin

RRAM Device Performances under Capacitive-Enhanced Current Programming Scheme

Keywords: resistive random access memory, RRAM, OxRAM, HfO2, operation scheme, current driver, current overshoot, switching variability

INTRODUCTION

RESISTIVE random access memories (RRAM) are emerging as a promising candidate for non-volatile memory (NVM) applications. Various memory suppliers are now bringing it up to the embedded market [START_REF] Peters | Reliability of 28nm embedded RRAM for consumer and industrial products[END_REF], [START_REF] Ito | ReRAM Technologies for Embedded Memory and Further Applications[END_REF]. Oxide based RRAM are made of a stack of metal / oxide / metal and are generally called OxRAM. The cells can be switched from a high-resistive state (HRS) to a low-resistive state (LRS) by electrical means; we call it the set operation. The reset operation makes the cells go back from LRS to HRS. Most of the set schemes are based on voltage drivers. During set operation, the abrupt reduction of the resistance induces a current surge. To prevent damaging the cell, voltage drivers incorporates current limiting elements. To enhance diverse metrics, refinements of this general scheme have been studied such as write termination [START_REF] Liu | 4.7 A 65nm ReRAM-enabled nonvolatile processor with 6× reduction in restore time and 4× higher clock frequency using adaptive data retention and self-write-termination nonvolatile logic[END_REF]- [START_REF] Chen | A 16Mb dualmode ReRAM macro with sub-14ns computing-in-memory and memory functions enabled by self-write termination scheme[END_REF], delayed write termination [START_REF] Yang | 24.2 A 14nm-FinFET 1Mb Embedded 1T1R RRAM with a 0.022µm2 Cell Size Using Self-Adaptive Delayed Termination and Multi-Cell Reference[END_REF], [START_REF] Chou | An N40 256K×44 embedded RRAM macro with SL-precharge SA and low-voltage current limiter to improve read and write performance[END_REF], and voltage ramps [START_REF] Sassine | Novel Computing Method for Short Programming Time and Low Energy Consumption in HfO2 Based RRAM Arrays[END_REF], [START_REF] Feng | Improvement of State Stability in Multi-Level Resistive Random-Access Memory (RRAM) Array for Neuromorphic Computing[END_REF].

The conductance in LRS is an important parameter as it is linked to metrics such as the read margin or the retention. While RRAM writing is field activated [START_REF] Ielmini | Evidence for Voltage-Driven Set/Reset Processes in Bipolar Switching RRAM[END_REF], the conductance in LRS has been reported to be strongly linked to the maximum current passing through the cell [START_REF] Ielmini | Resistive switching memories based on metal oxides: mechanisms, reliability and scaling[END_REF]. In light of this observation, it seems attractive to control the current rather than the voltage.

Current control writing has mainly been reported on in DC / quasi-static studies [START_REF] Li | Self-compliance SET switching and multilevel TaOx resistive memory by current-sweep operation[END_REF]- [START_REF] Lian | Improved Resistive Switching Uniformity in Cu/HfO2/Pt Devices by Using Current Sweeping Mode[END_REF]. To the best of the authors' knowledge, only a single publication [START_REF] Jain | 2 A 3.6Mb 10.1Mb/mm2 Embedded Non-Volatile ReRAM Macro in 22nm FinFET Technology with Adaptive Forming/Set/Reset Schemes Yielding Down to 0.5V with Sensing Time of 5ns at 0.7V[END_REF] wrote about using a current driver for writing RRAM at time scales of the order of the microsecond. However, its reliance on a voltage limiter subdues its potential benefits. The advantages and tradeoffs of current based programming against the voltage based approach were not clarified so far.

Another issue RRAM faces is high programming current. Since the main conducting mechanism in LRS is filamentary, the LRS conductance does not scale with the device [START_REF] Sandrini | OxRAM for embedded solutions on advanced node[END_REF]. This problem cannot be solved by programming the cells to a lower conductance since it is detrimental to its retention and to the read margin [START_REF] Chen | Improvement of data retention in HfO2/Hf 1T1R RRAM cell under low operating current[END_REF], [START_REF] Chen | Programming-conditions solutions towards suppression of retention tails of scaled oxide-based RRAM[END_REF].

In this letter, we explore a novel current driven writing scheme. To tackle the current scaling issue, we exploit the current spike caused by the inherent parasitic line capacitance in RRAM matrices. This enables this scheme to write cells to a high conductance with limited current amplitude. This helps relax issues with IR drop and pushes back the electromigration barrier, which are both concerns for scaling. It is also beneficial toward lowering the programming energy consumption.

EXPERIMENTAL

The measurements presented hereafter were obtained on the stack TiN / Ti / HfO2-5nm / TiN with a cell surface of 300 x 450 nm. We used a 130 nm CMOS technology with RRAM co-integrated in between M4 and M5 [START_REF] Sandrini | OxRAM for embedded solutions on advanced node[END_REF]. The RRAM are individually selected using pass-gates as we do not need current limitation in our proposed scheme.

RRAM Device Performances under Capacitive-Enhanced Current Programming Scheme

Thomas Bauvent, Paola Trotti, Olivier Billoint, Jean-François Nodin, Yasser Moursy, Gabriel Molas, Gaël Pillonnet, Univ. Grenoble Alpes, CEA, Leti, France Fig. 1c presents the operating principle of our writing scheme. As explained in the next section, with a current driver, the parasitic line capacitance Cp plays an important part in defining the dynamics of the system. To have a realistic capacitance value Cp with respect to Mb matrices, the current source must be co-integrated with the RRAM bank. With postlayout extraction, we estimated that Cp was approximately 600 fF in our integrated circuit. With scaling, and at constant bank size, this capacitance will tend to lower. Its value has no effect on the current spike amplitude, but is directly proportional to its energy [START_REF] Shrestha | Analysis and Control of RRAM Overshoot Current[END_REF]. This leads us to believe that there might be some need to adapt the scheme for more advanced technology nodes.

CAPACITIVE-ENHANCED CURRENT DRIVER

When a voltage driver is used, the parasitic line capacitance Cp does not play a major role. Other parasitic capacitances play a role in current spiking, such as the one on the node linking the cell to its current limiter [START_REF] Shrestha | Analysis and Control of RRAM Overshoot Current[END_REF]. However, when using a current driver, Cp plays an integral part in the writing process. That is why we call the driver in this scheme Capacitive-Enhanced Current Driver (CECD).

Fig. 2 presents a comparison between these two driving methods on a set operation. With voltage driving, when the set event occurs, the output current of the driver quickly rises to the compliance. As a result, the current limiter makes the voltage on the RRAM drop, but the current through the device stays at the compliance level until the voltage pulse ends.

In this work, with CECD, both the current passing through and the voltage on the RRAM rise at a limited rate which is proportional to Iprog Ú Cp where Iprog is the output current of the driver. When the set occurs, the voltage drops and the current passing through the RRAM spikes as Cp discharges through the RRAM. As long as the switching dynamics of the RRAM is faster than the RC time constant of the parallel combination of Cp and the impedance of the cell (≥ 5 kΩ x 600 fF = 3 ns in this work), the amplitude of the spike does not depend on the value of the capacitance [START_REF] Shrestha | Analysis and Control of RRAM Overshoot Current[END_REF]. The capacitance value only affects the total extra energy received. Since the conductance of the LRS is mainly defined by the maximum peak current passing though the RRAM device [START_REF] Ielmini | Resistive switching memories based on metal oxides: mechanisms, reliability and scaling[END_REF], we can expect the conductance in LRS to be well controlled, independently of the precise value of Cp. In the case of memory matrices, this should make this scheme robust to cell placement.

Thanks to the capacitive current spike, the maximum current passing through the RRAM is higher than Iprog. As such, with CECD, the output current of the driver can be significantly lower than the one in a voltage driver, for an equivalent LRS conductance. This reduces the constraints due to IR drop in memory matrices, pushes back the electromigration barrier, and reduces the total energy required per set.

RESULTS AND DISCUSSION

Fig. 3a shows the distributions of read currents obtained on a device for different Iprog. It shows that the distribution is mostly independent from the current pulse amplitude as long as it is above a given threshold, which in our experiment was found to be below 23.5 µA. In voltage driven schemes, a minimum voltage is required for writing the cells. Here, this limitation translates to a minimum current through the relation linking the voltage, current and impedance of the cell. To obtain the same conductance level with a typical voltage driver, a compliance current of at least 200 µA would be required. The irregularities near the high current tails seem to be linked to LRS instabilities that occur stochastically during cycling. This does not affect our conclusion on the independence of the LRS on Iprog, which is very contrasting compared to the importance of the compliance current in voltage schemes.

Fig. 3b shows the distributions of read currents obtained on a device for different pulse duration from 1 µs to 10 µs. Again, the read current distributions are found to be independent of the pulse duration, which is to be expected given that this scheme is naturally self-terminated. Indeed, by comparing the measured read current at 0.2 V with Iprog, we can infer that, after the set occurred, the voltage on the RRAM drops to a value close to 0.2 V. Therefore, from the point of view of the RRAM cell, the excitation roughly ends as soon as the spike from the set event wears off. That explains why the pulse duration has no measurable impact. Furthermore, the low remnant electric field could be beneficial to reduce relaxation as observed with delayed write termination [START_REF] Yang | 24.2 A 14nm-FinFET 1Mb Embedded 1T1R RRAM with a 0.022µm2 Cell Size Using Self-Adaptive Delayed Termination and Multi-Cell Reference[END_REF], [START_REF] Chou | An N40 256K×44 embedded RRAM macro with SL-precharge SA and low-voltage current limiter to improve read and write performance[END_REF].

Such as with the effect of the current pulse amplitude, we expect to hit a threshold in pulse duration under which the set operation would no longer be consistent. We successfully demonstrate that CECD works with pulses as short as 1 µs.

Since the energy consumption is proportional to both the pulse duration and current amplitude, shorter pulses would be beneficial not only for writing speed, but also for energy efficient programming. Fig. 3c shows a comparison of energy consumption between our scheme and [START_REF] Jain | 2 A 3.6Mb 10.1Mb/mm2 Embedded Non-Volatile ReRAM Macro in 22nm FinFET Technology with Adaptive Forming/Set/Reset Schemes Yielding Down to 0.5V with Sensing Time of 5ns at 0.7V[END_REF], [START_REF] Chou | An N40 256K×44 embedded RRAM macro with SL-precharge SA and low-voltage current limiter to improve read and write performance[END_REF], [START_REF] Sassine | Novel Computing Method for Short Programming Time and Low Energy Consumption in HfO2 Based RRAM Arrays[END_REF]. This establishes our scheme as an energy efficient writing technique that does not require any extra write termination circuitry.

Fig. 4a shows the scatter plot of the measured resistance under 0.2 V against the cycle index, in a 200 000 cycles experiment. After a period of seasoning of 2 000 cycles, both distributions remained stable throughout the remainder of the 200 000 cycles. Fig. 4b shows the distribution excluding seasoning. The observed raw bit failure rate is 5x10 -5 . For comparison, a distribution obtained with a conventional voltage scheme on the same technology is shown [START_REF] Giraud | Benefits of Design Assist Techniques on Performances and Reliability of a RRAM Macro[END_REF]. We can see that its LRS distribution is tighter, but achieves about the same raw bit failure rate. The low LRS/HRS ratio of less than one order of magnitude may seem insufficient for high capacity memories. However, by using circuit design strategies, a fully functional 128 kbit array with zero failures up to 1 M cycles was demonstrated on the same technology [START_REF] Giraud | Benefits of Design Assist Techniques on Performances and Reliability of a RRAM Macro[END_REF].

To further enhance the reliability, a write verify (WV) algorithm is employed. After each programming operation, if the resulting conductance is not satisfactory, the memory is cycled. This process is repeated until the conductance of the memory is within predefined bounds. Fig. 4c shows the histograms of the number of extra cycles required to meet the targeted read margin (RM). As expected, the larger the RM target, the more cycles are required. However, even for a target of 25 µA, only 4 extra cycles are required, at most, to meet our criteria on the 200 000 cycles of this experiment. In more than 99 % of the cases no extra cycle is necessary, and the maximum number of cycles is only necessary in rare events, less than 10 ppm. By comparison, the same write algorithm with the same rejection rate would only achieve a ~15 µA read margin with the conventional voltage scheme. For instance, a 25 µA read margin would hardly be possible with the voltage scheme. The lower slope of the distribution of LRS in the case of CECD, which could be seen as an inconvenience at first, actually makes WV more efficient. The WV algorithm requires extra cycles and as such consumes energy and ages the cell prematurely. However, since extra cycles are rare, its effect on global energy consumption and aging is marginal.

The low write energy exhibited by this scheme raises a question about retention. This work being the first on this new current scheme concept, we focused on its key and unique features. Therefore, we did not specifically address retention. We can report that we did not experience any issue with shortterm retention during this study. We plan to fully address retention in future work.

Given our previous observation that the LRS is independent of the value of Iprog above a threshold, the driver in this scheme need not be any more complicated than a fixed current source. As such, the surface required per driver is much lower than any voltage control implementation.

CONCLUSION

We presented a current driver based scheme for set operation in RRAM. The capacitive-enhanced current driver (CECD) utilizes the inherent parasitic line capacitance to set RRAM cells to a high conductance with limited current coming from the driver. We showed that the distribution of LRS is mostly independent of the pulse amplitude, duration and line capacitance. The scheme is used to program a cell for 200 000 cycles with a raw bit failure rate of 5 x 10 -5 . With write verification, it manages a read margin of 25 µA with only 4 extra cycles at worst. Since in more than 99 % of the cases the read margin is attained without requiring any extra cycle, the mean cost in time and energy of the write verification is marginal. The proposed scheme lowers the power consumption since it allows 10 x reduction in the current of the driver at a given LRS conductance compared to conventional drivers. This is a step forward to relax the constraints due to IR drop and electromigration and to promote the RRAM integration.

Fig. 1 .

 1 Fig. 1. (a) Photography of our test circuit. (b) Cross-section of the RRAM co-integrated in the BEoL (c) The writing scheme studied in this letter. The capacitance Cp models the parasitic line capacitance on BL in memory matrices.

Fig. 3 .

 3 Fig. 3. (a) Cycle-to-cycle distribution of read current for different write pulse current amplitude obtained with CECD. (b) Cycle-to-cycle distribution of read current for different write pulse duration. (c) Estimation of energy consumption of the CECD scheme compared to a current driver scheme ignoring the effect of the parasitic capacitance[START_REF] Jain | 2 A 3.6Mb 10.1Mb/mm2 Embedded Non-Volatile ReRAM Macro in 22nm FinFET Technology with Adaptive Forming/Set/Reset Schemes Yielding Down to 0.5V with Sensing Time of 5ns at 0.7V[END_REF], the extra energy invested in D-WT[START_REF] Chou | An N40 256K×44 embedded RRAM macro with SL-precharge SA and low-voltage current limiter to improve read and write performance[END_REF], and Ramp Voltage Stress (RVS)[START_REF] Sassine | Novel Computing Method for Short Programming Time and Low Energy Consumption in HfO2 Based RRAM Arrays[END_REF]. The difference of spread in the LRS distribution between (a) and (b) is due to the use of two different dies.

Fig. 2 .

 2 Fig. 2. A comparaison of the time evolution of the voltage and current in RRAM cells with a voltage driver on one side, and CECD on the other.

Fig. 4 .

 4 Fig. 4. (a) Scatter plot of the read resistance of 200 000 consecutive SET / READ / RESET / READ cycles with CECD. The bottom plot illustrates the effect of the WV with a RM of 25 µA. (b) Corresponding distribution (excluding the 2 000 seasoning cycles) with and without write verify (WV) and the distribution obtained in [19] with a conventional voltage scheme on the same technology. (c) Histograms of the number of extra cycles required by our WV scheme.

© 2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.


domain: phys spi

text:
François Richez

Rohit Jain

Prediction of the laminar-to-turbulent transition position on a helicopter rotor blade with cyclic pitch variation

This work aims at assessing different laminar-to-turbulent transition models with different CFD codes for unsteady cases of a helicopter rotor with cyclic pitch variation. ONERA evaluates with elsA code the semi-empirical transition criteria-based approach and the Langtry-Menter model, while US Army investigates two versions of the Langtry-Menter model available in the Helios/Overflow code. The experimental transition position measurements provided by DLR show a complex hysteretic behavior of the transition front motion. The simulations show that Langtry-Menter model can reproduce this hysteresis with a satisfactory agreement, if spatial accuracy is maximized using refined meshes and minimized numerical dissipation. The numerical results are then used to identify the transition mechanisms involved and the unsteady effects at the origin of hysteresis.

NOTATION

𝑎𝑎 ∞ freestream speed of sound, m/s b number of blades c blade chord, m 𝐶𝐶 𝑓𝑓 skin friction coefficient, 𝜏𝜏 𝑤𝑤 /[0.5𝜌𝜌 ∞ (𝑟𝑟Ω) 2 ] 𝐶𝐶 𝑇𝑇 /𝜎𝜎 rotor thrust coefficient, 𝑇𝑇/[𝜌𝜌 ∞ 𝑆𝑆𝜎𝜎(𝑅𝑅Ω) 2 ] 𝑀𝑀²𝐶𝐶 𝑛𝑛 section normal force in the local airfoil frame 𝑀𝑀 local Mach number, 𝑟𝑟Ω/𝑎𝑎 ∞ 𝑀𝑀 75 Mach number at 0.75𝑅𝑅, (0.75𝑅𝑅Ω)/𝑎𝑎 ∞ 𝑅𝑅 75 Reynolds number at 0.75𝑅𝑅, 𝜌𝜌 ∞ (0.75𝑅𝑅Ω)c/μ R rotor radius, m 𝑟𝑟 radial coordinate, m S rotor disk area, 𝜋𝜋𝑅𝑅 2 𝑇𝑇 𝑢𝑢 freestream turbulence level, % 𝑈𝑈 ∞ freestream axial velocity, m/s 𝑥𝑥 𝑡𝑡𝑡𝑡 transition position in the chordwise direction, m 𝜃𝜃 75 collective pitch angle at 𝑟𝑟/𝑅𝑅 = 0.75, ° 𝜃𝜃 1𝑠𝑠 cyclic longitudinal pitch angle, ° 𝜆𝜆 axial advance ratio, 𝑈𝑈 ∞ /(𝑅𝑅Ω) 𝜇𝜇 molecular viscosity, kg/m/s 𝜇𝜇 𝑡𝑡 eddy viscosity, kg/m/s 𝜌𝜌 ∞ freestream fluid density, kg/m³ 𝜎𝜎 rotor solidity, 𝑏𝑏𝑅𝑅𝑏𝑏/𝑆𝑆 𝜏𝜏 𝑤𝑤 wall shear stress, N/m² 𝜓𝜓 blade azimuth, ° Ω rotor angular velocity, rad/s

INTRODUCTION 1

Laminar-to-turbulent transition process can have a strong impact on helicopter rotor performance. In hovering flight, laminar boundary layer regions can be observed on a very large part of the lower side of the blade and on the leading-edge area of the upper side. Due to the reduced levels of skin friction in these laminar regions, the blade profile power is decreased which has a beneficial effect on the Figure of Merit. In forward flight, the blade sections see such large variations of both relative speed and local angle of attack that the transition position is expected to significantly vary with azimuth. However, the boundary layer transition is still often neglected in most of CFD simulations where the flow is usually assumed as fully turbulent. Consequently, the rotor torque is often overestimated with respect to experimental measurements. Turbulence transition models do exist but still suffer from a lack of validation for helicopter rotor applications, mainly because experimental measurements of boundary layer transition position on helicopter rotor blade are challenging and still few in the public literature.

Weiss et al. from the German Aerospace Lab (DLR) recently proposed high-quality transition position unsteady measurements based on Differential Infrared Thermography (DIT) technique on a small scale dynamically pitching rotor

Presented at the Vertical Flight Society's 79th Annual Forum & Technology Display, West Palm Beach, FL, USA, May 16-18, 2023. This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States. DISTRIBUTION STATEMENT A: approved for public release; distribution is unlimited. PR20230177. referred to as RTG hereafter (Ref. [START_REF] Weiss | Unsteady boundary-layer transition measurements and computations on a rotating blade under cyclic pitch conditions[END_REF]. CFD simulations were also performed in this article with the TAU code of DLR using URANS (Unsteady Reynolds-Averaged Navier-Stokes) modeling coupled with semi-empirical transition criteria approach. The experimental measurements show that the transition position as a function of the pitch angle exhibits a hysteresis loop, meaning that the transition position is different during the upstroke and the downstroke parts of the pitch cycle. The numerical simulation predicted well the transition position during downstroke but failed to accurately provide the transition during upstroke. One of the RTG test cases was later numerically investigated with transition modeling based on transport equations. Jain evaluated the Amplification Factor Transport (AFT) model and different versions of the Langtry-Menter model within Helios CFD code and showed very good prediction of the transition hysteresis loop (Ref. 2). Carnes and Coder also used the AFT model within OVERFLOW CFD code and paid a particular attention to the effect of freestream turbulence level (Ref. [START_REF] Carnes | Numerical Investigation of Unsteady Boundary Layer Transition on a Dynamically Pitching Rotor[END_REF]).

We propose here to push further the validation of the transition models by investigating two different test cases of the RTG database, and two different CFD solvers involving different transition models. After introducing the test cases and the numerical approaches, a numerical sensitivity analysis is performed in order to define the best numerical strategies. Then the numerical simulations are assessed by comparison with experimental data. In the last section, a discussion is conducted in order to clarify the transition mechanisms and the causes of hysteresis.

TEST CASES

The RTG rotor, defined in Table 1, is a four bladed rigid rotor with a radius R = 650 mm and a chord length c = 72 mm (Figure 1.a). The blades have a linear twist of -9.33°/R and are equipped with the DSA-9A airfoil (Figure 1.b). The rotor is operated with a small axial flow to prevent recirculation flow. The flow parameters of the two test cases considered here are presented in Table 2. Test case V gives a Mach number at r/R = 0.75 of 0.21 and cyclic pitch amplitude of 5.9° around a mean value of 9°. For test case II, the Mach number at 75% of radius is reduced to 0.11 while the cyclic pitch amplitude is decreased to 2.9° around a mean value of 10.1°. For both cases, the azimuthal variations of the transition positions on the upper side of the blade are measured by means of DIT technique as described in Ref. 1 (Figure 1.c). These experimental data will be used as reference for the validation of the CFD simulations.

NUMERICAL APPROACHES CFD solvers

Two different numerical approaches have been applied by the two partners of this joint research project.

US Army DEVCOM AvMC used the Helios software (Ref. 4). In Helios, Overflow (Ref. [START_REF] Nichols | OVERFOW User's Manual, Version 2.2[END_REF]) is involved for the blade body-fitted grids. Overflow used a Roe's 3rd-order upwind scheme using MUSCL (Monotonic Upstreamcentered Scheme for Conservation Laws) reconstruction for convective flux for both flow and turbulence. The 3rd-order flux used a 4th-order flux difference with a 3rd-order artificial dissipation term. Viscous fluxes were discretized using a second-order central discretization. The implicit time marching used the 2nd-order BDF (Backward Difference Formula). For temporal convergence, the Improved SSOR (Successive Symmetric Over-relaxation) scheme with 80 Newton iterations was employed. A time step size of 0.25° was used.

For the studies involving the effect of dissipation, the central difference scheme with different levels of dissipation was used for inviscid flux, as described later. For temporal discretization, the second-order diagonalized Beam Warming pentadiagonal scheme is used along with 80 dualtime stepping sub-iterations.

The off-body solver, SAMCART, is a structured, adaptive high-order Cartesian solver. The inviscid fluxes are discretized with a fifth-order central spatial scheme and a third-order explicit Runge-Kutta integration scheme is used for time advancement. The Helios simulations uses the Delayed-Detached Eddy Simulation (DDES) based on k-ω SST turbulence model coupled with Langtry-Menter (LM) transition model (Ref. 6). The initial LM model has been modified to guaranty Galilean invariance (Ref. 2) and to include crossflow transition correlation (Ref. 7). The Helios simulations based on Langtry-Menter with Galilean invariance is referred to as Helios-SST-LM-G and the simulations that also includes crossflow correlation are referred to Helios-SST-LM-G-CF in the following.

ONERA applied the elsA solver (Ref. 8) which discretizes the Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations with finite volume approach. The temporal derivative discretization involves a second-order implicit Gear scheme, with an azimuthal time step of 0.25° and 80 sub-iterations. All the simulations are performed with the k-ω SST turbulence model (Ref. 9). Two different turbulence transition models are evaluated. The first one applies semi-empirical transition criteria (TC) approach (Refs. 10 and 11) to detect laminar regions where the eddy viscosity is set to zero. The semi-empirical transition criteria are not computed locally but deduced from the chordwise history of the boundary layer, along the wall grid lines. Once the transition is detected by one of the criteria, an intermittency function is defined on the blade wall surface. Then, following the wall normal direction, the intermittency is used as a weight in the calculation of the eddy viscosity of the turbulence model, in order to set to zero the eddy viscosity in the laminar regions. This approach is thus not based on additional transport equations. The details of this TC model can be found in Ref. 11. The second approach evaluated with elsA is the Langtry-Menter (LM) transition model which is, for its part, based on transport equations of two new variables as described in Ref. [START_REF] Langtry | Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes[END_REF]. For TC simulations, the convective flux is discretized with Jameson's cell-centered second-order scheme using second-and fourth-order artificial viscosity coefficients (Ref. 12). For LM simulations, the convective flux is discretized with a specific version of the second-order AUSM+(P) scheme (Ref. [START_REF] Liou | 30 th Computational Fluid Dynamics, chap. AUSM schemes and extensions for low mach and multiphase[END_REF]) that permits to reduce the numerical dissipation for low Mach number flows (Ref. [START_REF] Mary | Large Eddy Simulation of Flow around an airfoil near stall[END_REF]. The elsA simulations performed with LM model and TC approach will be referred to as elsA-SST-LM and elsA-SST-TC hereafter.

Grid parameters

An overview of the grids used for the different solvers is shown in Figure 2. All the simulations are based on overlapping grid strategy. The parameters of the blade grid and the background grid are detailed in Table 3 andTable 4 respectively.

With Helios (Figure 2.a), the blade grids are composed of a main grid of 14.1 million nodes, completed with overlapping tip and root caps, leading to a total number of 22 million nodes per blade. The center body is also included in the simulation with a body-fitted grid of 11.3 million nodes. The background Cartesian grid is composed of several blocks with several levels of resolution ensuring the finest grid resolution of 5% of chord around the rotor. It has 78 million points.

For elsA simulations, the rotor model is more simplified. First, the center body is not included as shown in Figure 2.b. Secondly, only the profiled part of the blade is meshed while the blade root is not considered. Two blade grids were built in order to investigate the sensitivity of the solution with respect to the grid resolution. The first grid (Figure 3.a), composed of 5.8 million nodes, has the same resolution as the grids classically used for rotor simulations at ONERA. It is referred to as "coarse grid" hereafter. The second grid (Figure 3.b) is significantly refined leading to a total number of 15.5 million points per blade. The background Cartesian grids are composed of 19 million nodes for the coarse grid and 48 million nodes for the fine grid (Table 4).

All the parameters of the blade grids used by the partners are compared in detail in Table 3. One can remark that the grid resolution in the chordwise direction is very similar between Helios grid and elsA fine grid, while elsA coarse grid is approximately coarsened by a factor two. In the spanwise direction, Helios grid is slightly finer than the two elsA grids. In the wall normal direction, the resolutions are similar between all grids, especially the wall normal spacing that has been chosen in order to ensure Δy + ≤ 1.

Extraction of the transition position

Estimating the transition location from a numerical simulation based on Langtry-Menter model is not straightforward. The value of the eddy viscosity ratio or the intermittency is imposed to zero at wall. Therefore, an extraction of the flow field on the blade surface cannot be useful to detect the transition. Thus, computing the transition position requires extracting a turbulent quantity at a certain distance from the wall that is to be defined by the user. Since the three CFD codes used here does not offer the same post-processing functionalities, each partner has defined his own strategy to estimate the transition position.

On ONERA side, the extraction of the transition position is done as follows: the flow field on the blade grids is stored at each time step. Then, a surface is extracted at a constant grid index in the wall normal direction (k). The value of the k-index is chosen to correspond approximately to the logarithmic region of the boundary layer. Then, at several blade sections, the transition position is defined as the first grid point in the chordwise direction where the eddy viscosity ratio gets larger than a user-defined threshold. Since the k-index could be out of the logarithmic region of the boundary layer, the process is repeated for five different values of k and the final transition position is defined as the minimum value of all the results. The transition location in Helios/Overflow computations is determined using a special post-processing of the boundary layer flow field, as follows. For each azimuth and each radial station, boundary layer over the blade section is examined for growth in turbulence. The chord location on the section where the turbulence level in the boundary layer has grown by an order of magnitude with respect to the inflow boundary is flagged as the transition location. Since structured grids were used for blades, this post-processing method could be applied in a straight-forward manner as the grid lines are aligned along the surface tangent and normal directions.

NUMERICAL SENSITIVITY

In order to define the best numerical strategy, the partners first investigated the sensitivity of the solution with respect to numerical parameters. In this preliminary activity, only Case V of the database is considered (Table 2).

ONERA investigated the sensitivity with respect to the time integration with elsA solver. For that purpose, the number of sub-iterations of the Newton process used to solve iteratively the non-linear system at each time step is increased from 20 to 40 and finally 80, while the azimuthal time step is kept to a constant value of Δ𝜓𝜓 = 0.25°. All the simulations are performed with LM model on the coarse grid defined in Table 3 and shown in Figure 3.a. The transition position 𝑥𝑥 𝑡𝑡𝑡𝑡 is defined as the chordwise position where the eddy viscosity ratio gets higher than 𝜇𝜇 𝑡𝑡 /𝜇𝜇=1. The azimuthal evolutions of 𝑥𝑥 𝑡𝑡𝑡𝑡 thus defined are compared at the section 𝑟𝑟/𝑅𝑅 = 0.75 in Figure 4.a. As expected, all the simulations predict that the transition 𝑥𝑥 𝑡𝑡𝑡𝑡 /𝑏𝑏 is located around 80% of chord in the first part of the cycle (0°≤ 𝜓𝜓 ≤ 180°) where the pitch angle is the highest, while the transition moves to the leading edge region in the second half of the cycle (180°≤ 𝜓𝜓 ≤ 360°) where the pitch angle reaches its maximum value. However, it can be seen that the azimuthal angle where the transition moves upstream and downstream is very sensitive to the number of sub-iterations used in the time integration. In particular, the upstream transition motion occurs at 𝜓𝜓 = 280° with 20 sub-iterations instead of 𝜓𝜓 = 240° when 80 sub-iterations are used. This high sensitivity with respect to the time integration is probably due to the fact that the transition motions seem to occur in a very fast manner. In light of this observation, it will be chosen to fix the number of sub-iterations to a minimum value of 80 in the rest of this work. It can also be remarked that the transition position shows very high oscillations when it is located upstream, i.e. for 250°≤ 𝜓𝜓 ≤ 360°. This will be discussed in the next paragraph.

The influence of the grid resolution is also investigated for the LM model and still with the elsA code. Two simulations are performed for Test Case V, using the coarse and the fine grids shown in Figure 3 and defined in Table 3. The variations of the transition position at 𝑟𝑟/𝑅𝑅 = 0.75 are compared in Figure 4.b. The grid refinement perfectly eliminates the large oscillations of the transition position 𝑥𝑥 𝑡𝑡𝑡𝑡 /𝑏𝑏 observed for 250°≤ 𝜓𝜓 ≤ 360° and previously mentioned. The reason lies in the nature of the transition process. Indeed, as indicated by the sign of the skin friction coefficient at 𝑟𝑟/𝑅𝑅 = 0.75 and 𝜓𝜓 = 270° (Figure 5) the transition occurs in a Laminar Separation Bubble (LSB) that extends along 15% of chord. The fine grid allows to discretize the LSB with 100 grid points in the chordwise direction while the coarse grid only provides less than 20 points inside the LSB. The poor resolution of the coarse grid does not permit a sufficient production of turbulence kinetic energy to ensure a stable reattachment point of the LSB (Figure 6). Thus, the very high resolution of the fine grid as defined in Table 2 seems mandatory to properly capture the transition process for this specific test case. The grid refinement also slightly shifts downstream the position of the transition for 0°≤ 𝜓𝜓 ≤ 180° but this effect is less significant as the one observed in the second half of the cycle. As a conclusion, the fine grid will be considered in the following, since it provides a better capture of the LSB. The influence of the numerical dissipation was also assessed using Helios/Overflow. For this study a central difference scheme was used for inviscid flux discretization as then accuracy and dissipation levels can be independently varied. Table 5 shows the various settings. The first row, DISS1, is the baseline case, which used a Roe's 3rd-order upwind scheme. Other cases, DISS2, DISS3, and DISS4 used a central difference scheme in the order of increased numerical dissipation -DISS2 used a fifth-order scheme with a small artificial dissipation coefficient of 0.01; DISS3 had the coefficient increased to 0.04; and DISS4 was most dissipative with second-order central with a large dissipation coefficient. The effect on transition location prediction is illustrated in Figure 7. DISS1 and DISS2 have minimal dissipation and the predictions agree closely. However, as the dissipation is increased the transition onset is delayed, in general. The Helios/Overflow results presented hereafter used the DISS1 settings.

As previously mentioned, extracting the transition location 𝑥𝑥 𝑡𝑡𝑡𝑡 from a numerical simulation based on LM model is not straightforward. In the 𝑥𝑥 𝑡𝑡𝑡𝑡 calculation method presented above, the threshold value of 𝜇𝜇 𝑡𝑡 /𝜇𝜇 used to define the transition has to be arbitrary set by the user. In order to evaluate the sensitivity of 𝑥𝑥 𝑡𝑡𝑡𝑡 with respect to this threshold value, the transition position obtained with the criteria of 𝜇𝜇 𝑡𝑡 /𝜇𝜇 = 1 and 0.1 were compared in elsA computations. In both cases, transition positions are computed from the same flow field provided by elsA code for the Test Case V. The two threshold values are significantly lower than the level reached in turbulent boundary layer. Thus, we can expect that 𝑥𝑥 𝑡𝑡𝑡𝑡 will be weakly sensitive to these values. However, the results shown in Figure 8 indicate that the effect is not negligible. The threshold 𝜇𝜇 𝑡𝑡 /𝜇𝜇 = 0.1 provides a transition location that is located upstream by 2% of chord for 240°≤ 𝜓𝜓 ≤ 360° and 5% of chord for 60°≤ 𝜓𝜓 ≤ 230°. These results show that the chordwise evolution of turbulence production, and consequently of eddy viscosity, is rather slow compared to more usual high Reynolds number flows. From the following, it is decided to keep a threshold of 𝜇𝜇 𝑡𝑡 /𝜇𝜇 = 0.1 since it provides a more continuous evolution of 𝑥𝑥 𝑡𝑡𝑡𝑡 especially in the azimuthal portion 0°≤ 𝜓𝜓 ≤ 60°.

VALIDATION

Following the best practices defined in the previous section, US Army and ONERA have performed the simulations for both Test Case V and Test Case II.

Test Case V

The rotor maps of the transition positions for Test Case V are shown in Figure 9. It is reminded in the experimental rotor maps (Figure 9.a) that the minimum pitch angle is reached at ψ=90° and the maximum pitch angle at ψ = 270°. The rotor map of the experimental measurements also reveals that the transition position is not symmetrical with respect to lateral axis. In particular the upstream motion of the transition occurs around 210° of azimuth and the downstream motion occurs around 30° of azimuth. This means that the transition position evolves differently between the upstroke and downstroke parts of the pitch cycle.

The rotor maps of transition position provided by Helios simulations (Figure 9.b and Figure 9.c) show a good agreement with respect to the experiment. The addition of the crossflow correlation into the original LM transition model slightly improves the results, especially at the inboard sections around ψ = 210° and at the outboard sections around ψ = 30°. Despite these slight differences, one can conclude that both models provide the good trend, especially the dissymmetry of x tr between upstroke and downstroke parts of the pitch cycle.

The elsA simulation based on the LM model is shown in Figure 9.d. Unlike LM model in Helios, no correction has been added to original model in elsA code in order to ensure the Galilean-invariance. However, the elsA-SST-LM transition map (Figure 9.d) looks rather similar to the one provided by Helios (Figure 9.b). Once again, the asymmetry of the rotor map is quite well predicted, although the upstream transition motion during upstroke and the downstream transition motion during downstroke are slightly delayed.

The elsA simulation based on the TC approach provides the rotor map shown in Figure 9.e. The agreement with experiment is not as good as with LM model. The downstream motion of transition around 30° of azimuth is quite well predicted. However, the strong delay of the upstream motion of x tr during the upstroke is not captured. Finally the TC approach gives a symmetric transition map, which indicates that this transition model does not predict significant differences of x tr between upstroke and downstroke.

For a more qualitative assessment of the models, the azimuthal evolution of the transition position at r/R = 0.75 is shown in Figure 10 for both codes. Figure 10.a confirms that Helios SST-LM-G provides a very good agreement with respect to experiments. Including crossflow transition correlation into the model (SST-LM-G-CF) permits to avoid a spurious evolution of x tr observed for 20°≤ ψ ≤ 40° but also slightly degrades the agreement with experiment for 90°≤ ψ ≤ 180°. The elsA simulation with LM model shown in Figure 10.b provides a prediction similar to Helios-SST-LM-G. In particular, the spurious evolution of x tr between 20° and 40° of azimuth is again observed. The transition position is however a little more upstream with elsA for all azimuths. Another significant difference between the LM model's results of elsA and Helios lies on the upstream motion of transition that occurs between 180° and 270° of azimuth. While Helios prediction is perfectly in phase with the experimental measurements, elsA shows a noticeable delay of the upstream motion of x tr . The reason of the discrepancy between Helios-SST-LM-G and elsA-SST-LM can be manifold. As seen in the previous section, the solution seems very sensitive to the numerical parameters. Thus, the differences of numerical methods between Helios and elsA could explain the discrepancy. The corrections introduced in Helios to allow for Galilean invariance, that are not used in elsA, could also be the reason of this slight disagreement between codes. Looking now at the elsA results based on TC approach in Figure 10.b, one can observe a good agreement of the transition position during downstroke, i.e. from ψ = 270° to ψ = 90°. However, during upstroke, the TC model fails to predict the strong delay of the upstream transition motion.

Plotting the transition position as a function of the pitch angle allows to highlight the dissymmetric between upstroke and downstroke in the form of a hysteresis loop. As shown in Figure 11, the hysteresis observed with DIT measurements is large. For instance, at 𝜃𝜃 75 = 9°, the transition is located at 30% of chord during downstroke and 80% of chord during upstroke. LM model with both Helios and elsA code permits to reproduce this large hysteresis, athough the agreement is less good with elsA in the second part of downstroke. TC approach predicts a small hysteresis, narrow around the experimental measurements during downstroke.

Test Case II

As shown in Table 2, Test Case II differs from Test Case V by a lower amplitude of the cyclic pitch (2.9° instead of 5.9°), a slightly higher collective pitch (10.1° instead of 9°) and a lower rotor angular velocity which leads to a decrease of the Reynolds number at 0. 75𝑅𝑅 from 3.16 × 10 5 to 1.72 × 10 5 .

The Helios and elsA results for Test Case II are shown in Figure 12 in the form of transition maps. The experimental measurements reveal again a dissymmetry of 𝑥𝑥 𝑡𝑡𝑡𝑡 between upstroke and downstroke. Helios simulations based on LM model (Figure 12.b and Figure 12.c) show a satisfactory prediction of the transition motions although the agreement is not as good as for Test Case V. The crossflow correlation (Helios-SST-LM-G-CF) improves the prediction especially during upstroke, for 0°≤ 𝜓𝜓 ≤ 30°, similarly to what was observed for Test Case V. The elsA simulation with LM model (Figure 12.d) gives the good trend but shows the same shortcoming as already observed for Test Case V. Indeed, the upstream motion of transition for 180°≤ 𝜓𝜓 ≤ 270° is too delayed. Furthermore, the downstream motion of transition that occurs between 0° and 30° of azimuth shows a disrupted evolution. With TC approach (Figure 12.e), the transition seems located to upstream at all azimuths, and more significantly for 30°≤ 𝜓𝜓 ≤ 210°.

A detailed look to 𝑟𝑟/𝑅𝑅 = 0.75 (Figure 13) confirms the observations previously mentioned. Helios provides a good agreement of the transition between 0° and 220° of azimuth. The upstream motion of 𝑥𝑥 𝑡𝑡𝑡𝑡 is also well-captured, but the transition is located too upstream for 230°≤ 𝜓𝜓 ≤ 360°. With elsA code, the LM model gives an unexpected "bump" of 𝑥𝑥 𝑡𝑡𝑡𝑡 between 45° and 90° of azimuth, similarly to what was observed for Test Case V, but this spurious behavior is here more pronounced. Again, the upstream motion of transition (180°≤ 𝜓𝜓 ≤ 270°) is delayed and too abrupt in comparison with the experimental measurements. At the beginning of the downstroke (270°≤ 𝜓𝜓 ≤ 360°) the transition position is very similar to Helios results but located too upstream in comparison with the experiment. The TC approach shows a less satisfactory prediction of 𝑥𝑥 𝑡𝑡𝑡𝑡 . For 0°≤ 𝜓𝜓 ≤ 180°, the transition is located too upstream. It reaches 60% of chord instead of 85% in the experiment. On the second half of the cycle, the transition position shows some oscillations probably due to a lack of turbulence production in the LSB.

Plotting the transition position with respect to the pitch angle 𝜃𝜃 75 reveals that Test Case II leads to a hysteretic behavior of the transition (Figure 14) in the same way as Test Case V. As observed for Test Case V, the LM model simulations predict this large hysteresis while the TC approach underpredicts its amplitude.

DISCUSSION

We propose in this section to analyze the origin of the hysteresis observed in the evolution of the transition location as a function of the pitch angle. Since the phenomenon looks similar between the two test cases, we will focus the discussion on Test Case V only. A first possible reason of the transition hysteresis could be a hysteresis in lift. Indeed, if lift, and thus pressure distribution, is different during upstroke and downstroke, the transition position will consequently differ. Figure 15 shows the evolution of the section normal force at 𝑟𝑟/𝑅𝑅 = 0.75 as a function of the pitch angle 𝜃𝜃 75 predicted by numerical simulations. One can see that both Helios-SST-LM-G-CF and elsA-SST-TC provide the same normal force evolution as the one provided with a fully turbulent simulation (referred-to as "elsA-SST-FT" in the figure). Furthermore, the hysteresis in normal force is very weak, which suggests that the hysteresis of transition could not be due (or at least only due) to a hysteresis in lift. This confirms what is suggested by Weiss et al. in Ref. [START_REF] Weiss | Unsteady boundary-layer transition measurements and computations on a rotating blade under cyclic pitch conditions[END_REF]. The semi-empirical transition criteria used for elsA-SST-TC simulation are based on quasi-steady assumptions. At each time step, the transition position is deduced from the pressure distribution and the properties of the boundary layer at this current time step only. Thus, the TC approach of elsA code is able to predict hysteretic behavior of transition only if this one is caused by hysteresis in lift. Since the hysteresis in lift is here small, it is not surprising to observe that elsA-SST-TC simulation predicts a small transition hysteresis. Furthermore, the discrepancy between measurements and elsA-SST-TC results finally give more credence to the idea that the transition hysteresis observed in the experiments is not caused by lift hysteresis.

Additional elsA-SST-TC simulations have been performed without cyclic pitch and with a collective pitch angle of 6°, 9° and 12°. The transition positions obtained with these three "non-cyclic" simulations are added to the Test Case V results in Figure 16, in order to highlight how the unsteady case differs from a quasi-steady condition. The figure indicates that the transition position of Test Case V is very close to the non-cyclic results during upstroke, while it strongly diverges from them during upstroke. This tends to prove that the transition motion during downstroke follows a quasi-steady state, and that the hysteresis is mainly due to unsteady effects during upstroke.

In order to identify these unsteady effects, we propose in Figure 17 and Figure 18 to follow the azimuthal evolution of the transition position together with the skin friction distribution 𝐶𝐶 𝑓𝑓 at 𝑟𝑟/𝑅𝑅 = 0.75. Figure 17 shows Helios-SST-LM-G-CF and elsA-SST-LM results that are assumed as representative of the hysteresis observed experimentally. Figure 18 shows the elsA results using TC approach that are considered as the quasi-steady solution.

At the beginning of downstroke, (𝜓𝜓 = 270° and 𝜓𝜓 = 315°), the transition occurs in a LSB located around 20% of chord, as indicated by the negative value of 𝐶𝐶 𝑓𝑓 . Both LM model and TC approach provide almost the same transition positions and skin friction distribution.

On the second part of downstroke, from 𝜓𝜓 = 0° to 𝜓𝜓 = 90°, all the transition models predict that the transition progressively moves downstream to finally reaches 80% of chord. However, the differences observed in the skin friction evolution between Figure 17 and Figure 18 points out that the transition process is different between the two models. On one hand, LM model (Figure 17) predicts that transition is still induced by a LSB from 0° to 90° of azimuth and that the LSB moves downstream together with the transition point. On the other hand, the transition predicted by TC approach (Figure 18) occurs at the limit of flow separation at 𝜓𝜓 = 0°. Then, we observe an attached boundary layer transition that moves downstream until 𝜓𝜓 = 90°.

During the first part of upstroke (from 90° to 210° of azimuth), the quasi-steady assumptions (Figure 18) shows that the attached flow transition moves progressively upwards, preventing any laminar separation in the leadingedge region. With LM model (Figure 18), in contrast, the transition is stuck in the trailing edge region, which permits another laminar separation to emerge around 20% of chord at 𝜓𝜓 = 210°. Between 210° and 215°, LM model show two flow separation regions, a first one close to the trailing edge which corresponds to the vanishing LSB responsible of transition during the previous azimuths, and a second one around 20% of chord which corresponds to the arising LSB responsible of the transition at the current time. It is noted that prediction of the coexistence of the two LSBs in Helios simulation is limited to a small range of azimuth (210°≤ 𝜓𝜓 ≤ 215°), while in elsA-SST-LM simulations it is stretched up to 𝜓𝜓 = 232°. However, in both cases, the coexistence of the two separation regions is observed.

In the end of upstroke, from 240° to 270° of azimuth, the trailing edge flow separation disappeared and all the simulations whether with LM model or TC approach gives a LSB between 10 and 30% of chord.

CONCLUSIONS

Different turbulence transition models have been evaluated for unsteady conditions of a small rotor in lowspeed axial flow with cyclic blade pitch variations. Two CFD solvers have been used: Helios on US Army side and elsA on ONERA side. US Army has evaluated two versions of the Langtry-Menter model, the standard version based on Langtry correlation and a second one including crossflow transition correlation. Since the Langtry-Menter model is not Galilean invariant, some corrections were introduced in Overflow to allow the invariance for both versions.

ONERA has evaluated methods for transition modeling: the Langtry-Menter model (without Galilean invariance) and a semi-empirical transition criteria approach.

A numerical sensitivity analysis has first been performed and the conclusions are as follows:

• The Langtry-Menter model is very sensitive to the grid resolution, especially in the chordwise direction. In the region where laminar separation bubble appears, a too coarse grid will lead to a lack of turbulence production and the appearance of non-physical oscillations of the reattachment point. The grid resolution required for the RTG test cases are significantly higher than the one typically used for full-scale helicopter rotor flow simulations.

• The Langtry-Menter model is also very sensitive to the numerical dissipation and the accuracy of the time integration. These parameters mainly impact the lag in azimuthal phase of the upstream motion of the transition front.

• The Transition-Criteria approach does not require such refined grid and small time-step to get a converged solution (although it is not shown in the paper).

The capabilities of the simulations to predict the unsteady motions of the transition front have then been assessed and the analysis of this assessment has led to the following conclusions:

• The experimental measurements reveal a strong hysteretic behavior of the transition motion, meaning that the transition positions are different during upstroke and downstroke parts of the pitch cycle.

• LM model in Helios code provides a very good agreement with respect to the experiment, especially concerning the hysteresis. The crossflow transition correlation slightly improves the prediction.

• LM model in elsA code provides a satisfactory agreement with experiment but not as good as Helios.

The slight discrepancies between the two codes could be due to numerical methods or Galilean invariance corrections.

• The Transition-Criteria approach provides the good estimation of the upstream and downstream maximum position of transition for Test Case V but not for Test Case II. In both cases, this approach fails to predict the large hysteresis of the transition.

Finally, an analysis has been conducted in order to understand the transition mechanisms and the origins of this hysteresis. It has led to the following conclusions:

• The hysteresis in transition is probably not caused by hysteresis in lift.

• The transition position during downstroke follows a quasi-steady state, while unsteady effects appear during upstroke.

• Since the transition-criteria approach is based on quasisteady assumptions, it can only reproduce the downstroke part of the cycle.

• According to LM model, the transition is induced by a laminar separation bubble almost all along the pitch cycle.

• The lag phase of the transition position during the upstroke leads to the coexistence of two laminar separation bubbles during a short range of azimuth.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Overview of the RTG rotor test model (from Ref. 1).

Figure 5 .Figure 6 .

 56 Figure 5. Skin friction coefficient at 𝒓𝒓/𝑹𝑹 = 𝟎𝟎. 𝟕𝟕𝟕𝟕 and 𝝍𝝍 = 𝟐𝟐𝟕𝟕𝟎𝟎° (LM model, elsA, fine grid).

Figure 7 .

 7 Figure 7. Sensitivity of the Langtry-Menter model with respect to numerical dissipation (Helios simulations).

Figure 8 . 13 Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .

 81391011121314 Figure 8. Sensitivity of 𝒙𝒙 𝒕𝒕𝒓𝒓 calculation with respect to the eddy viscosity ratio threshold 𝝁𝝁 𝒕𝒕 /𝝁𝝁.

Figure 15 .

 15 Figure 15. Section normal force 𝑴𝑴 𝟐𝟐 𝑪𝑪 𝒏𝒏 as a function of pitch angle at 𝒓𝒓/𝑹𝑹 = 𝟎𝟎. 𝟕𝟕𝟕𝟕 (Test Case V) including the results for non-cyclic conditions with elsA-SST-TC (red triangle symbols).

Figure 16 .

 16 Figure 16. Transition position as a function of pitch angle at 𝒓𝒓/𝑹𝑹 = 𝟎𝟎. 𝟕𝟕𝟕𝟕 (Test Case V) including the results for noncyclic conditions with elsA-SST-TC (red triangle symbols).

Figure 17 .

 17 Figure 17. Skin friction coefficient 𝑪𝑪 𝒇𝒇 at 𝒓𝒓/𝑹𝑹 = 𝟎𝟎. 𝟕𝟕𝟕𝟕 for Test Case V: elsA-SST-LM (green line), Helios-SST-LM-G-CF (blue line). Letter "T" indicates the transition position for elsA (green) and Helios (blue).

Figure 18 .

 18 Figure 18. Skin friction coefficient 𝑪𝑪 𝒇𝒇 at 𝒓𝒓/𝑹𝑹 = 𝟎𝟎. 𝟕𝟕𝟕𝟕 for Test Case V with elsA-SST-TC. Letter "T" indicates the transition position.

Table 1 . Rotor blade parameters.

 1

	Radius, 𝑅𝑅	650 mm
	Chord length, 𝑏𝑏	72 mm
	Tip sweep	Parabolic (0.91 ≤ 𝑟𝑟/𝑅𝑅 ≤ 1)
	Twist rate	-9.33°/𝑅𝑅
	Airfoil	DSA-9A
	Number of blades, 𝑏𝑏	4
	Solidity, 𝜎𝜎	0.141

Table 2 . Flow parameters of the two test cases of the RTG database.

 2

	Case V	Case II

Table 3 . Blade grid parameters.

 3

	Helios	elsA coarse grid	elsA fine grid

Table 4 . Background Cartesian grid parameters.

 4

	Helios	elsA coarse grid elsA fine grid

Table 5 . Numerical dissipation parameters evaluated with Helios.

 5

		Order	IRHS	Dissipation	Number of
				Coefficient	Sub-iterations
	DISS 1	Third	Roe Upwind	-	80
	DISS 2	Fifth	Central	0.01	80
	DISS 3	Fifth	Central	0.04	40
	DISS 4	Second	Central	0.04	40

ACKNOWLEDGMENTS

The authors would like to acknowledge the U.S./France Project Agreement on Rotary-Wing Aeromechanics and Human Factors Integration Research, a longstanding cooperation between the U.S. Army, ONERA and DGA.

The authors would like to acknowledge the DLR for sharing the experimental measurements of the RTG rotor test that were used for the validation of the numerical simulations.


A Gassenq
email: alban.gassenq@univ-lyon1.fr

H-S Nguyen

E Cleyet-Merle

S Cueff

A Pereira

Diffraction grating enhanced photoluminescence from etching-free erbium thin films

Micro-structuration by etching is commonly used in integrated optics, adding complex and costly processing steps which can also potentially damage the devices performances due to the degradation of the etched sidewalls. For diffraction grating fabrication, different strategies have been developed to avoid etching, like layer deposition on structured surface or grating deposition on top of active layers. However, etching remains one of the best process for making high aspect ratio diffraction grating. In this work, we have developed fully structured diffraction gratings (i.e. like fully etched gratings) using lift-off based processing performed in pulsed laser deposited layers since the combination of both techniques is of great interest for making microstructures without etching. We have first studied the influence of the lithography doses in the lift-off process, showing that (1) micrometric spatial resolution can be achieved and (2) sidewall angle can be controlled from 50 to 150° in 0.5 µm thick layers. Using such optimizations, we have then fabricated Er-doped Y2O3 uniaxial diffraction gratings with different periods ranging from 3 to 8 µm. The fabricated devices exhibit emission and reflectivity properties as a function of the collection angle in good agreement with the modeling with a maximum luminescence enhancement of x15 compared to an unstructured layer at a wavelength of 1.54 µm. This work highlights thus the lift-off based processing combined with pulsed laser deposition as a promising technique for etch-free practical applications such as luminescence enhancement in Er doped layers.

INTRODUCTION

Material etching is a standard process in photonics for the fabrication of a large variety of micro-devices like waveguides, lasers, photodetectors gratings, and so on. It consists of the ablation of a surface trough a micro-structured mask (e.g photoresist). Such a processing step remains complicated, expensive and potentially detrimental for practical applications since it can also induce extra-losses due to the degradation of materials along etched sidewalls [1,2] as well as undesired roughness. For instance, grating devices are widely spread in photonics for light coupling [3], filters [4], mirror definitions [5] and so one. Furthermore, luminescence enhancement by diffraction grating is also very interesting to increase the spontaneous emission, to control the light direction or to increase the pump absorption [6][7][8][9][10][11]. For that purpose, different strategies have been developed to avoid etching in the emissive layer like using diffraction elements below [6][7][8]12] or above the active layer [9,10], but these methods cannot be used to fully structure the emitting layers (like the fully etched gratings compared to the shallow etched gratings). For rare earth doped materials, etching is also commonly used in photonics [13][14][15][16][17]. Other methods have also been studied for integration but with millimetric resolution [18,19], without structuration [10,20,21], using ion beam [22][23][24], or by layer deposition on pre-structured substrates [25,26]. Indeed, rare earth emitters have recently growing as a very interesting emerging platform for quantum photonics [27,28] due to the lanthanide high coherences [29] allowing to realize many quantum functions like sources [16], spin readout [17,24], memories [22], microwave transductors [23] and so on. Therefore, developing diffraction gratings with rare earth emitters without etching is of great interest to control the light emission of such layers for different applications like laser integration [20,30] with grating, signal processing [31], quantum optics [27,28], and so on.

In this work, we have fabricated Er 3+ rare-earth-doped gratings by Pulsed Laser Deposition (PLD) and lift-off processing without etching. Lift-off processing is usual in microelectronics but is limited to low thermal budget layers [32,33] and can also induce high sidewall roughness [34]. Since PLD deposition can provide high quality films at low temperature deposition due to the high kinetic energy of the atoms during deposition [35], it is thus a compatible technique with the lift-off processing to fabricate thicker micro-devices without etching [36]. However, low roughness and high-resolution patterning in thick layer are relatively challenging with lift-off processing. For that purpose, we have first investigated the lithography dose influence in the shape of the microstructure to determine the minimum resolution size. We have then fabricated diffraction gratings with different periods and measured their emission and reflection properties as a function of the wavevector. Structured layers exhibit a controlled angular-dependence in emission, in good agreement with theory and a strong luminescence enhancement at certain angles. This work thus shows that the lift-off processing combined with PLD is a relevant approach for etch-free fabrication of practical photonic devices like diffraction gratings as future building blocks in many applications in integrated optics like rare-earth ions based quantum nanophotonics. [27]

GRATTING

A. Description

Our reference design consists of 0.6 µm thick Y2O3:Er 3+ layer on SiO2 with 1.87 [37] and 1.44 [38] refractive index at 1.54 µm wavelength, respectively. This structure corresponds to a single-mode waveguide slab in the infrared spectral window 1.3-1.8 µm. Most of Y2O3:Er 3+ emission radiates into this guided mode and is trapped within the slab due to total internal reflection, while only ~8% percent of the emitted light may escape to the air. To enhance light extraction and absorption, a periodic corrugation of period p is introduced (see Figure 1-a). This modulation diffracts the initial guided mode and induces a band folding effect: each diffraction order m implements a wavevector-shift of Km=2mπ/p to the dispersion of the initial guided mode. As a consequence, for certain values of diffraction order m, the folded guided mode is brought above the light-line and becomes guided mode resonances (GMR). The GMRs are out-coupling channels for Y2O3:Er 3+ guided emission and the angular pattern of emitted light is dictated by the dispersion of the GMRs. At first approximation, when only the bandfolding mechanism is taken into account, the dispersion of GMRs can be analytically calculated. The gray shading region indicates the numerical aperture in our experimental setup (NA=0.42) studied later in this work. For a small period of p=2µm, we "only" have eight GMRs corresponding to diffraction orders m=±1, ±2, ±3, ±4. When limited within a numerical aperture NA=0.42 and focusing in the vicinity of erbium emission wavelength at 1.54µm, only two GMRs of m= ±2 exists for 2 µm period which are responsible of the extraction of erbium emission to free space. However, for big periods, many more GMRs are present, and the light extraction mechanism involves multi channels. Similar arguments can be applied for the enhancement of absorption mechanisms in the visible range (λp=0.532 µm for the 4 I15/2→ 2 H11/2 transition of Er 3+ ). In that case, the waveguide is already multimode transverse in the visible range, thus the number of GMRs in the visible range is much more important (bandfolding for different guided modes, with different diffraction orders). Therefore, the mode positions have been summarized as a function of the period at 0° on the Figure 1-b around the Er 3+ extracted emission (λ=1.54 µm for the 4 I13/2 → 4 I15/2 transition) and absorption (λp=0.532 µm for the 4 I15/2→ 2 H13/2 transition) which will be experimentally studied later in this work. We clearly see that the lower period absorbs at lower diffraction mode (with different transverse modes k=1, 2 or 3) and that the emission is mono-mode (with only one transverse mode k=1).

B. Modeling

FABRICATION

A. Processing

Figure 2 presents the used process flow for the grating fabrication. Micro-strips were first patterned with AZ5214 photoresist previously spin coated on top of 2 µm thick SiO2 layer on Si with a laser writer using GDS files. Different exposure doses were investigated, whereas the development time was kept constant at 60 s. A 500 nm-thick Y2O3:Er 3+ (1 at. %) layer was then deposited on top of the structured substrate with a KrF excimer laser (λ=248 nm, t=17 ns, Coherent Compex Pro) operating at 5 Hz. The laser beam was focused on the target over 2 mm 2 , and the laser energy and the Y2O3:Er 3+ substrate distance were kept constant at 73 mJ and 6 cm, respectively. Substrate rotation was used to obtain homogeneous layer thickness at a growth speed of 3.5 nm/min. The film was grown at room temperature in an oxygen-gas atmosphere (10-3 mbar). Lift-off was then performed in ultrasonic wet bench with acetone. After the lift-off, samples were annealed for 24 h at 650 °C in order to activate the doping and to improve the structural quality of the film [36]. c) with SEM (Scanning Electron Microscopy) images. SEM images were taken at 45° angle tilt with a TESCAN SEM at 5 kV with a secondary electron detector. Low doses were obtained with negative photoresist and high doses with positive photoresist using AZ5214E[Eu] reversible photoresist from micro-chemicals company, in order to control the slope angle of the photoresist [39]. By comparing the imaging before and after lift-off, we clearly see that the PLD layer follows the lithography shape. Figure 3-b plots the slope θ angle of the microstructure and the pattern width difference Δa (see inset of the Figure 3-b) as a function of the lithography dose. This width difference corresponds to the absolute difference of the pattern size compared to the theoretical width (a in the Figure 1-a), which was used in the GDS file and sent to the laser lithography. For negative photoresist, the pattern increases with the dose since it remains on the sample after exposure (i.e. negative mode). For positive photoresist, the dimension of the developed pattern decreases with the dose since the exposed region is removed from the sample after exposure (i.e. positive mode). Regarding the shape of the devices, the obtained microstructures have a sidewall angle ranging from 50° to 150°. For positive photoresist, the sidewall angle is below 90° with higher roughness coming from the breaking of the layer during the lift-off. For negative photoresist, low roughness is observed since the deposited material is not in contact with the photoresist, especially for the lowest dose. These behaviors come from the reversible nature of the used photoresist AZ5214. Therefore, the exposure dose is a primordial parameter which needs to be controlled as a function of the targeted shape. The Figure 3-c presents two examples of extrema for the minimum and maximum angles we can obtain on a disk. In that case, we observe here "bowl" (for low dose with negative photoresist) or "pancake" (for high dose with positive photoresist) shapes which confirms that we can obtain atypical shapes and very low roughness (for low dose) with this technique. Indeed, as it can be seen in the SEM images after liftoff for low dose in the Figure 2-a, the sidewall roughness is in the same order of magnitude than the roughness of the unstructured part of the layer.

FIG. 3: Influence of the lithography dose on the shape of the microstructure with a) SEM tilted imaging before and after lift-off; b) sloped sidewall angle and width difference (compared to the theoretical width of the used GDS file); c) Examples of the disk micro-structuration at very low and very high doses, resulting in "bowl" and "pancake" shapes, respectively.

GRATTING CHARACTERIZATION

A. Measurements

Different gratings have been fabricated with 0.5 µm thick Y2O3:Er 3+ on 2 µm SiO2 on Si substrate. Figure 4-a presents SEM images for 3, 4.5, 6 and 8 µm periods. We clearly see that high spatial resolution is obtained with sub-micrometric trenches. Photoluminescence and reflectivity measurements were performed on such gratings. Figure 4-b presents the micro-reflectivity measurements performed with an infrared lamp focused over 10 µm² on the microstructures. Light was collected through a x50 objective of numerical aperture NA=0.42 and directed to a spectrometer with an IR InGaAs camera used as output detector with 10 ms integration time. In order to resolve the angular distribution, the farfield image is captured in the Fourier space by projecting the backfocal plane of the objective on the entrance slits of the spectrograph. The corrugation direction in the sample is aligned along the spectrometer slit, so that single-shot image taken by the camera provides directly an angle-resolved reflectivity map [40]. The GMRs are revealed in the reflectivity as Fano resonances in the measured spectrum. The associated modeling (similar to the one discussed in Figure 1) have been added as red dashed lines on the measurements with the adjusted period pth indicated on top of the Figure 4-b. We found a very good agreement between the measured and adjusted periods starting from 4 GMRs for the 3 µm period up to 10 for the largest period, as expected. Micro-PL measurements were also performed on the same microstructures with a 532 nm wavelength CW laser and the same measurement setup (Figure 4-c). First of all, we observe the expected erbium transition in Y2O3 poly-crystalline matrix with a maximum emission at 1.54 µm wavelength [41]. Secondly, we observe that the photoluminescence intensity is highly increased at angular positions driven by the GMR dispersions. For instance, the emission is enhanced at two angles (±10°) for the 3 µm period, which correspond to the diffracted modes observed by the reflectivity measurement (Figure 4).

B. Photoluminescence enhancement

In order to evaluate quantitatively the luminescence enhancement, the gratings were compared to the unstructured part of the sample. The Figure 5-a presents the luminescence of the bulk layer with the same measurement setup. Since the emitted signal was lower, the integration time used was 10 times higher (i.e 100 ms). We also see a maximum of emission at 1.54 µm. However, no preferential orientation of the emission is observed since rare earth emission is isotropic in bulk materials. Such emission was compared in count/s with the structured layers at 1.54 µm as a function of the emitted angle (Figure 5-c) and at -10° as function of the wavelength (Figure 5-d). A maximum luminescence enhancement of 15 is reported for the 3 µm period compared to the bulk layer at the maximum of the detected signal (i.e. 1.54 µm and -10°). Such improvement results from the extraction enhancement evaluated to be a factor 2 (by comparing the intensity between 0 and ±10°) and the absorption enhancement evaluated here to be around 7.5 (15 divided by the extraction factor of 2). Indeed, the absorption enhancement of the pump laser is a key factor for luminescence enhancement in diffraction gratings [10]. As it was shown in Figure 1-b, at 0.54 µm wavelength the coupled in-plane modes have an order of around 15 for the 3 µm period while it is closer to 40 for the 8 µm period. According to the diffraction grating theory, the diffraction will be more efficient with lower period since the diffraction occurs with lower mode giving thus higher absorption. Our maximal enhancement factor (i.e. x15) could thus be further increased using lower periods. Furthermore, the luminescence enhancement is also related to the extraction factor. As indicated by the comparison between the reflection measurement and the photoluminescence (Figure 4-b and c), for the 3 µm period sample, we clearly see that the crossing of both Fano resonances does not fit the maximum of the erbium emission at 1.54 µm. Therefore, the previous enhancement value of 15 at λ=1.54 µm could also be further improved by carefully adjusting the Fano resonances at the maximum of the erbium emission using period optimization. Finally, far field measurement was performed on the 6 µm period sample using a spectral filter and a Fourier lens on the beam path.

As expected three main modes are observed only along the grating diffraction (Figure 5-b), which confirms that we can control the emission along different directions with such a design.

FIG. 1 :

 1 FIG. 1: Diffraction grating modeling with (a) used parameters, (b) mode position as a function of the period around λEr=1.54 µm (top) and λp=0.532 µm (bottom) wavelength and (c) calculated dispersion of GMRs.

Figure 1 -

 1 Figure 1-c presents calculated GMR dispersions for the 1.3-1.8 µm wavelength ranges, with three different periods 2, 5 and 8 µm. The red dashed line corresponds to the wavelength at which the erbium emission is maximal (λEr=1.54 µm for the 4 I13/2 → 4 I15/2 transition of Er 3+ ). The gray shading region indicates the numerical aperture in our experimental setup (NA=0.42) studied later in this work. For a small period of p=2µm, we "only" have eight GMRs corresponding to

FIG. 2 :Figure 3 -

 23 FIG. 2: Process flow for the grating fabrication using a) lithography of photoresist, b) Y2O3:Er 3+ PLD, c) lift-off and d) annealingB. Dose influenceFigure3-a presents the influence of the lithography dose before and after the lift-off processing (corresponding to the process step in Figure2-b and 2-c) with SEM (Scanning Electron Microscopy) images. SEM images were taken at 45° angle tilt with a TESCAN SEM at 5 kV with a secondary electron detector. Low doses were obtained with negative photoresist and high doses with positive photoresist using AZ5214E[Eu] reversible photoresist from micro-chemicals company, in order to control the slope angle of the photoresist[39]. By comparing the imaging before and after lift-off, we clearly see that the PLD layer

Fig. 4 .

 4 Fig. 4. Diffraction gratings characterization for different periods with a) SEM tilted image; b) reflectivity and c) photoluminescence measurements as a function of the wavelength and the angle.

Fig. 5 .

 5 Fig. 5. Diffraction gratings with a) photoluminescence of the unstructured layers as a function of the wavelength and the angle; b) far field photoluminescence measurement in the Fourier plane; Comparison of the emission intensity for the different periods compared to the unstructured layers c) as a function of the angle and d) wavelengthCONCLUSIONIn summary, we have investigated Y2O3:Er 3+ diffraction gratings made by PLD and lift-off processing. The processing conditions have been first optimized showing that the lithography dose is a key parameter to control the sidewall angle and to reach high-resolution patterning. Diffraction gratings have then been fabricated down to a period of 3 µm. Reflectivity and photoluminescence measurements of such devices show a very good agreement with modeling, with a maximum luminescence enhancement of 15 compared to unstructured layer. Therefore, this work highlights the PLD based lift-off processing to fabricate photonics devices. It opens thus perspectives to the realization of new microstructures for integrated optics without etching, for future applications in free-space communications and integrated quantum optics.

Acknowledgments. The authors would like to thank the NanoLyon platform for the cleanroom facilities. Salahedine Toubi, Ylyes Betka for their M1 projects and the ILM for funding.

Disclosures. The authors declare no conflicts of interest.


text:
Rocco Moretti

Hyeonsoo Yeo

Franc ¸ois Richez

Rotor Loads Prediction on UH-60A Flight Test using Loose Fluid/Structure Coupling

This work exploits the high-quality UH-60A flight test campaign for comparison with concurrent simulations at the U.S. Army and ONERA. Rotor airloads and structural loads are predicted by coupling Computational Fluid Dynamics (CFD) with rotorcraft comprehensive analysis (CA). A high-speed test point is examined. Comparisons are shown for airloads, structural loads and rotor controls. Agreement with test data is improved when the flexibility of the pitch link is considered. Results show good agreement between the two partners predictions and with test data.

INTRODUCTION

In the last decades, numerical simulation has taken its place in the study of aerodynamic loads on aircraft and helicopters. It is essential to reveal and resolve possible design deficiencies/inconsistencies before wind tunnel or flight tests are performed.

In the case of helicopters, accurate rotor loads prediction is a really challenging task because of the interaction of several complex phenomena: unsteady low-speed and transonic aerodynamics of flexible blades undergoing large displacements and rotations. This work validates rotor loads prediction with existing tools at ONERA and U.S. Army on the UH-60A flight test data (Ref. [START_REF] Bousman | UH-60A Airloads Catalog[END_REF]. Similar comparison has already been carried out on ONERA 7A rotor (Refs. [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor Using Loose Fluid/Structure Coupling[END_REF][START_REF] Yeo | High-Fidelity Structural Loads Analysis of the ONERA 7A Rotor[END_REF]. The UH-60A rotor has already been extensively studied by the U.S. Army (Refs. [START_REF] Yeo | Rotor Structural Loads Analysis Using Coupled Computational Fluid Dynamics/Computational Structural Dynamics[END_REF][START_REF] Potsdam | Rotor Airloads Prediction Using Loose Aerodynamic/Structural Coupling[END_REF].

The state-of-the-art in the rotor loads prediction is coupling the rotorcraft comprehensive analysis (CA) with the computational fluid dynamics (CFD) (Ref. 5). Indeed, this allows to improve fidelity of aerodynamics by replacing the CA loworder aerodynamic model with the most accurate flow solver for this kind of computations.

The novelty of this work stands on four points: [START_REF] Bousman | UH-60A Airloads Catalog[END_REF] The cooperation between the U.S. Army Combat Capabilities Presented at the Vertical Flight Society's 79th Annual Forum & Technology Display, West Palm Beach, FL, USA, May 16-18, 2023. This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States. DISTRIBU-TION STATEMENT A: approved for public release; distribution is unlimited. The test data used in the present study were obtained during the NASA/Army UH-60A Airloads Program conducted from August 1993 to February 1994 (Figure 1). This seminal flighttest program provided a detailed airloads and blade structural loads database for a full-scale rotor covering a broad flight envelope. The UH-60A main rotor is a four-bladed, articulated rotor. The blade construction uses a titanium spar with a fiberglass outer contour and two airfoils, the SC1095 and SC1094R8 (Figure 2).

Development Command

The present study investigates a high-speed flight condition: counter c8534, µ = 0.37, C T /σ = 0.081, with advancing blade negative lift.

COMPUTATIONAL METHODOLOGIES

Rotor loads prediction is accomplished by coupling a rotor comprehensive analysis with a CFD code solving the Unsteady Reynolds-averaged NavierStokes equations (URANS). The CFD airloads replace, through a delta method, the air- U.S. Army used RCAS (Ref. 6) for comprehensive analysis and Helios (Ref. 7) for CFD. RCAS is a comprehensive multidisciplinary computer software system for predicting rotorcraft aerodynamics, performance, stability and control, aeroelastic stability, loads, and vibration. RCAS is capable of modeling a wide range of complex rotorcraft configurations operating in hover, forward flight, and maneuvering conditions. RCAS has been used extensively for correlation of performance and loads measurements of the UH-60A in various flight conditions. Helios is the rotary-wing product of the U.S. Army and CREATE Air Vehicles program sponsored by the U.S. Department of Defense High Performance Computing Modernization Office. Helios uses an innovative dualmesh paradigm that employs unstructured meshes in the near body close to the surface to capture the wall-bounded viscous effects and structured Cartesian grids in the off body to resolve the wake through a combination of higher-order algorithms and adaptive mesh refinement (AMR). The unstructured solver NSU3D is used for near-body and the structured solver SAMARC is used for the Cartesian off-body grid system. An overset procedure facilitates the data exchange and also enables the relative motion between the two meshes. The parallel domain connectivity solver PUNDIT automatically handles the data exchange between the two meshes.

ONERA used HOST (Ref. 8) for comprehensive analysis and elsA (Ref. 9) for CFD. HOST (Helicopter Overall Simulation Tool) is a rotorcraft comprehensive analysis developed by Airbus Helicopters. HOST modeling of blade dynamics is multibodylike. The blade is represented by a series of rigid elements connected by virtual hinges. Inertial properties (e.g., mass and moment of inertia) are assigned to each rigid element, and elastic properties (e.g., flap bending, chord bending, and torsion stiffness) are assigned to each hinge. A modal reduction approach is used to reduce the number of degrees of freedom from a large system of equations. The aerodynamics of HOST is based on a lifting-line approach based on airfoil tables combined with a wake model. In this effort, among the several wake models available, a prescribed wake helical geometry was used. The elsA CFD code, developed at ONERA, solves the unsteady Reynolds-averaged Navier Stokes equations for both background Cartesian grids and blade curvilinear grids. An overset procedure facilitates the data exchange and also enables the relative motion between the two meshes. The near-body grids of the blades are rotated and deformed following the blade motion and trim provided, through the loose coupling, by the rotorcraft comprehensive analysis HOST.

In the current work, ONERA and U.S. Army used similar grid meshes. Helios used unstructured near-body grids, whereas elsA used structured near-body grids, see Figure 3. The four rotor blade grids in Helios have 15.4 million nodes and 36.5 million cells, whereas in elsA the cell count is 16 million. The finest off-body spacing in Helios is 5% chord with a fixed refinement region surrounding the rotor plane. The off-body grid contains 146 million unblanked grid points on 8 levels. In elsA, the finest off-body spacing is 10% chord. Table 1 summarizes the main parameters characterizing the rotor and the studied flight point.

RESULTS

The results chapter is organised in two parts: the first part shows how the structural modeling of the UH-60A blade in HOST was modified in order to model the pitch link stiffness Once the structural models are shown to be equivalent, the second part compares the airloads, structural loads and control angles.

Campbell diagrams

Structural dynamics modeling of the UH-60A rotor blade is briefly described. First are compared the Campbell diagrams when considering a rigid pitch link in HOST and the nominal pitch link stiffness in RCAS, see Figure 5. It can be seen that there is a good agreement for the lower modes, while significant differences arise in the torsion mode.

The flexibility of the pitch link is then introduced in HOST by lowering the torsional stiffness of the blade at the radius corresponding to the attachment of the pitch link. This process was done iteratively, judging by the agreement between RCAS and HOST. Eventually, a torsional stiffness adjustment was obtained such that an equivalent pitch link behaviour was created in HOST, yielding a very good agreement with the RCAS Campbell diagram, see Figure 6.

As a conclusion, this chapter has shown the relevance of modeling properly the pitch link stiffness in the rotorcraft comprehensive analyses. The remainder of the results presented in this paper correspond to coupled CFD/CA analyses and include this modeling.

Airloads

In this work calculated airloads were obtained by integrating pressure all over the chord, and not only at the location of the pressure taps. This is expected to be acceptable given that up to 30 unsteady pressure transducers were installed at every instrumented section of the UH-60A blade (Ref. [START_REF] Biedron | An Examination of Unsteady Airloads on a UH-60A Rotor: Computation versus Measurement[END_REF]), yet any error in the blade pressures may change significantly the integrated pitching moments. Consequently, all plots of pitching moment have the mean removed. Pitching moment is marked by a pronnounced negative, nosedown, pitching moment on the advancing side, resulting from the highly compressible flow on the outboard part of the blade. Agreement between both partners calculations and test data is good, although less so for the most outboard section, at r/R=0.965.

Trim control angles

The rotor control angles used for trimming the simulation are the collective pitch angle, the lateral cyclic pitch angle and the longitudinal cyclic pitch angle. Rotor shaft angle of attack is prescribed at the experiment value, -7.31 deg. Figure 10 compares the calculated and measured oscillatory chord bending moment at r/R=0.113, 0.30, 0.50 and 0.70. Here the differences between the two coupled anlayses (elsA/HOST and Helios/RCAS) are more marked than for the flap bending moment, although both capture the experimental trends. Figure 13 compares the half peak-to-peak chord bending moment, showing for elsA/HOST a closer match to experiment. The significant underprediction of half peak-topeak chord bending moment was previously observed by both RCAS and CAMRAD II (Ref. 4).

Figure 11 compares the calculated and measured oscillatory torsion moment at r/R=0.30 and 0.70. Modelling a rigid or flexible pitch link in HOST has a visible influence on the predicted torsion moment, with the rigid pitch link exhibiting higher amplitude oscillations and fewer harmonic content. Figure 14 compares the half peak-to-peak torsion moment and shows that the equivalent pitch link model substantially improves the correlation with the test data and both partners underestimate the torsion moment amplitude in the outboard part of the blade.

Finally, Figure 15 compares the harmonic magnitude of calculated and measured structural loads. In this comparison, elsA/HOST analysis results with rigid pitch link are not included for clarity. This comparison is analysed next in the light of the following assumptions: harmonics 1 and 2 are more related to rotor trim, lag damper and pitch-lag coupling. harmonics of chord bending moment are harder to capture; elsA/HOST shows a quadratic trend, whereas Helios/RCAS a more linear one. Measurements are somewhere in-between. However, the elsA/HOST worsening agreement with experiment at the blade root suggests a difference in the way the lag damper is modelled between HOST and RCAS. The fourth harmonic of the chord bending moment remains ellusive to both calculations, whereas the third and fifth are better captured by Helios/RCAS. Overall, compared to chord bending moment, the harmonics of flap bending moment and torsion moment exhibit better agreement between calculations and with measurement.

The challenges in capturing chord bending moment were not observed in (Refs. 2, 3), where calculations were compared to the 7A rotor wind tunnel experiment. This may be related to the fact that rotor-drivetrain coupling may be smaller in the 7A rotor wind tunnel rig than in the UH-60A flight testing. It should be noted that the effects of drivetrain dynamics on the chord bending moment for the flight test and the effects of hub impedance using test stand NASTRAN model on the chord bending moment for the wind-tunnel test are also investigated for the UH-60A rotor. In general, the effects of drivetrain dynamics slightly improved both 4/rev and 5/rev harmonic magnitude correlation and the effects of hub impedance improved the correlation of both 3/rev and 4/rev harmonic components. Readers who are interested in those details are referred to (Ref. 11).

CONCLUSIONS

This work is an achievement of the cooperation between the U.S. Army and ONERA under the U.S./France Project Agreement on Rotary-Wing Aeromechanics and Human Factors Integration Research, a more than 50 year old cooperation between the United States of America and France in the field of rotorcraft.

Results support previous evidence that accounting for pitch link flexibility is an enabling factor for high-fidelity predictions of rotor loads.

The good agreement between the predictions of the two partners paves the way towards joint investigations of more challenging flight conditions and blade shapes, in order to ultimately improve the performance and loads prediction capability of high-performance, high-speed rotors.

  Aviation & Missile Center (DEV-COM AvMC) and the ONERA French Aerospace Research Lab for the analysis of the UH-60A rotor under the United States/France Project Agreement on Rotary-Wing Aeromechanics and Human Factors Integration Research, (2) Success of ONERA's coupled CFD/CA methodology on the prediction of the UH-60A rotor loads, (3) Systematic comparison of two high-fidelity CFD/CA coupled analyses, and (4) Demonstration of how simple assumptions on blade modeling in the CA tool can have an important impact on rotor loads prediction.

Figure 1 .

 1 Figure 1. UH-60A flight test picture.

Figure 2 .

 2 Figure 2. UH-60A rotor blade planforms showing locations of airfoil section, pressure taps and strain gauges. loads of the CA. The coupling exchanges are done on a per revolution (periodic) basis (loose coupling).

Figure 3 .

 3 Figure 3. Structured overset grids.

Figure 4

 4 shows the RCAS modeling of the rotor hub with blade hinges, pitch control, and the lag damper. The blade is composed of 13 nonlinear beam elements. Three coincident hinge elements are offset from the center of rotation. The three hinge elements allow for simultaneous flap, lag, and pitch rotations of the blade. Rigid bars and spring elements are used to represent the pitch control linkage.HOST assumes, by default, the pitch link to be rigid (i.e., infinite stiffness). An equivalent pitch link model was developed in this work by using torsionally soft beam element. More recently a pitch control elasticity model has been introduced in HOST by its main developers at Airbus Helicopters, but the authors did not have the time to use it for this paper. Detailed blade frequency comparisons are shown below. Campbell diagrams are calculated in vacuum conditions.

Figure 4 .

 4 Figure 4. RCAS modeling of the UH-60A rotor.

Figure 5 .

 5 Figure 5. Fan plot; rigid pitch link in HOST.

Figures 7 and 8 Figure 6 .

 86 Figures 7 and 8 compare the normal force and pitching moment, respectively, at three span-stations. The black line corresponds to measurement. Blue, to Helios/RCAS calcula-

Figure 9

 9 Figure 9 compares the calculated and measured oscillatory flap bending moment at r/R=0.113, 0.30, 0.50 and 0.70. The two coupled anlayses (elsA/HOST and Helios/RCAS) agree well with each other and with experiment. The equivalent pitch link model reduces amplitude at r/R=0.7 compared to the rigid pitch link model, and thus improves the correlation with test data. Figure 12 compares the half peak-to-peak flap bending moment. The rigid pitch link model significantly overpredicts the half peak-to-peak amplitude at the outboard locations. The equivalent pitch link model improves the correlation by decreasing the half peak-to-peak amplitude. Due to strong couplings among the high-frequency modes shown in Figure 6, pitch link stiffness has an important influence on flap bending moment. Overall, elsA/HOST tends to overpredict the amplitudes between r/R=0.40 and r/R=0.70, while the opposite is true for Helios/RCAS.

Figure 15 .

 15 Figure 15. c8534. Comparison of harmonic magnitude of calculated and measured structural loads.

Table 1 .

 1 Summary of UH-60A rotor flight test

	Characteristic	Symbol	Value
	Rotor Geometry
	Rotor radius	R	26.833 ft (8.179 m)
	Mean blade chord	c	1.722 ft (0.525 m)
	Number of blades	N b	4
	Rotor solidity	σ	0.0817
	Airfoils		SC1095,SC1094R8
	UH-60A test point C8534: High speed
	Air density	ρ inf	0.002082 slug/ f t 3
	Air density (metric)	ρ inf	1.0732 kg/m 3
	Air temperature	T inf	71.8 F (295.3 K)
	Rotational speed	Ω	258.1 rev/min
	Advance Ratio	µ	0.37
	Rotor thrust coefficient	C T /σ	0.081
	Sideslip angle	β	1.27 deg

and hence become comparable to RCAS. This was achieved by modifying the inboard torsional stiffness in HOST and comparing the Campbell diagrams produced by HOST and RCAS.

Table 2

 2

	com-

Army predicts an angle nearly 2 degrees smaller (in absolute value) than measurement, while ONERA underestimation is close to 3 degrees. Previously published work by the U.S. Army (Ref.

[START_REF] Yeo | Rotor Structural Loads Analysis Using Coupled Computational Fluid Dynamics/Computational Structural Dynamics[END_REF]

showed Helios/CAMRAD II to calculate very similar pitch control angles compared to Helios/RCAS for point c8534.

Table 2 .

 2 UH-60A c8534 control angles Collective Lateral cyclic Long. cyclic

		θ 0	θ 1c	θ 1s
	c8534	13.73	2.84	-11.05
	Helios/RCAS	13.61	1.81	-9.35
	elsA/HOST	13.87	3.38	-8.23
	Structural loads

  Harmonics higher than 2 are dominated by the blade elastic model and the coupling with the drivetrain. On the two first harmonics, the U.S. Army calculations (Helios/RCAS) feature better agreement with experiment, particularly for flap bending moment and torsion moment. The two first Figure 7. Sectional normal force coefficient C n M 2 .

		0.25		r/R=0.675 (elsA r/R=0.652)	0.25		r/R=0.865 (elsA r/R=0.860)	0.25		r/R=0.965 (elsA r/R=0.964)
		0.20					0.20					0.20
		0.15					0.15					0.15
		0.10					0.10					0.10
	n M 2	0.05					0.05					0.05
	C
		0.00					0.00					0.00
		0.05					0.05					0.05
		0.15 0.10	0	45	90 c8534 elsA/HOST 135 180 225 270 315 360 azimuth, deg elsA/HOST rigid pitch link Helios/RCAS	0.15 0.10	0	45	90	135 180 225 270 315 360 azimuth, deg	0.15 0.10	0	45	90	azimuth, deg 135 180 225 270 315 360
		0.010 0.015		r/R=0.675 (elsA r/R=0.652) elsA/HOST elsA/HOST rigid pitch link Helios/RCAS c8534	0.010 0.015		r/R=0.865 (elsA r/R=0.860)	0.010 0.015		r/R=0.965 (elsA r/R=0.964)
		0.005					0.005					0.005
	mean removed	0.005 0.000					0.005 0.000					0.005 0.000
	C m M 2	0.010					0.010					0.010
		0.015					0.015					0.015
		0.020					0.020					0.020
		0.025	0	45	90	135 180 225 270 315 360 azimuth, deg	0.025	0	45	90	135 180 225 270 315 360 azimuth, deg	0.025	0	45	90	azimuth, deg 135 180 225 270 315 360

Figure 8. Sectional pitching moment coefficient C m M 2 (mean removed).

ACKNOWLEDGMENTS

The authors would like to acknowledge the U.S./France Project Agreement on Rotary-Wing Aeromechanics and Human Factors Integration Research, a longstanding cooperation between the U.S. Army, ONERA and DGA.