V Jaunet

C Braud

Experiments on lift dynamics and feedback control of a wind turbine blade section

Keywords: Feedback control, Wind turbine, Lift

An experimental campaign is performed to the study of feedback lift control applied to a wind turbine blade. A 5-digit NACA profile whose trailing edge is rounded for circulation control purposes is used in combination with fluidic actuation. It is first shown in this article how the modified profile performs, in terms of aerodynamic forces, in both natural and manipulated cases. Then, the dynamics of controlled pressure (and thus lift) establishment is identified.

A discussion is given on the speed of lift dynamics in comparison with previous studies and on its ability to overcome lift perturbations due to external perturbations. Finally, a feedback lift control experiment is performed showing the feasibility of such control in a wind turbine environment.

Introduction

Wind turbines evolve in the atmospheric boundary layer at a height where strong wind velocity fluctuations in both amplitude and direction can be found [START_REF] Kaimal | Spectral characteristics of surface-layer turbulence[END_REF]. These fluctuations not only modify the actual forces acting on the blade but can also lead to dynamic stall [START_REF] Choudhry | An insight into the dynamic stall lift characteristics[END_REF][START_REF] Melius | Dynamic stall of an experimental wind turbine blade[END_REF], therefore creating undesirable mechanical Recently, Choudhry et al. [START_REF] Choudhry | Methods to control dynamic stall for wind turbine applications[END_REF] proposed to manipulate the flow using passive actutator to reduce the lift pertubation due to dynamic stall. Although the use of passive actuators revealed to be efficient in this purpose, it may be difficult to optimize their design for a variety of flow regimes. Moreover, the fact that they have to be permanently installed may cause efficiency penalty at standard functionning regime that might not be accepted by manufacturers and energy producers. Flaps and slats are also considered for this purpose [START_REF] Castaignet | Model predictive control of trailing edge flaps on a wind turbine blade[END_REF][START_REF] Castaignet | Full-scale test of trailing edge flaps on a vestas v27 wind turbine: active load reduction and system identification[END_REF][START_REF] Ortega Gmez | Load reduction potential on large wind turbines using active trailing edge flaps[END_REF][START_REF] Neuhaus | Mitigating loads by means of an active slat[END_REF] but additional moving parts, even very efficient, are often discarded by manufacturers due to the plausible additional maintenance costs. On the other hand, Active Flow Control (AFC) devices, such as systems using air injection to manipulate lift, has received a lot of attention in aeronautic applications to achieve a lighter and mechanically simpler high lift systems for transonic commercial transport aircraft [START_REF] Hartwich | Afc-enabled simplified high-lift system integration study[END_REF].

Since the early work of Englar and Williams [START_REF] Englar | Design of a circulation control stern plane for submarine applications[END_REF] and Englar [START_REF] Englar | Two-dimensional subsonic wind tunnel tests on a cambered 30 percent-thick circulation control airfoil[END_REF][START_REF] Englar | Experimental investigation of the high velocity coanda wall jet applied to bluff trailing edge circulation control airfoils[END_REF], circulation control (CC) using secondary injection and Coanda effect at the trailing edge is known to produce an efficient way of modifying the lift of a wing profile.

The lift gain produced by Coanda jet circulation control is obtained through two mechanisms. The main one is the displacement of the stagnation points at the leading edge and trailing edge, hence changing the circulation around the wing [START_REF] Englar | Experimental investigation of the high velocity coanda wall jet applied to bluff trailing edge circulation control airfoils[END_REF]. The second one is through suction by the very low wall pressure that the Coanda turning can produce on the Coanda surface [START_REF] Kadosch | Déviation des jets par adhérence à une paroi convexe[END_REF]. Circulation control using Coanda effect has been investigated for a great variety of field from helicopter blades to propellers, showing the great interest of the research community in such techniques [START_REF] Kweder | Applications of circulation control, yesterday and today[END_REF]. However, circulation control using Coanda has only recently been investigated in the field of wind turbine aerodynamics [START_REF] Dumitrache | Active and passive circulation control as enhancement techniques of wind turbines performance[END_REF][START_REF] Shires | Application of circulation controlled blades for vertical axis wind turbines[END_REF][START_REF] Baleriola | Circulation control on a rounded trailing-edge wind turbine airfoil using plasma actuators[END_REF][START_REF] Leroy | Comparison of flow modification induced by plasma and fluidic jet actuators dedicated to circulation control around wind turbine airfoils[END_REF][START_REF] Braud | Jet flow control at the blade scale to manipulate lift[END_REF]. Cagle and Jones [START_REF] Cagle | A wind tunnel model to explore unsteady circulation control for general aviation applications[END_REF] and Jones et al. [START_REF] Jones | An active flow circulation controlled flap concept for general aviation aircraft applications[END_REF] provided experimental feasability evidence of CC on wing profiles similar to the one used on wind turbines but only gave steady-state information on the system. To the author's knowledge, few studies aimed at revealing the dynamics of lift enhancement by means of circulation control systems while this is an essential Closed-loop control of the blade lift must be considered to dynamically alleviate the mechanical loads due to wind gusts. Closed-loop lift control has been of great interest especially for post-stall performance recovery of aircraft airfoils. Williams et al. [START_REF] Williams | Lift response of a stalled wing to pulsatile disturbances[END_REF][START_REF] Williams | Closed-loop control of a wing in an unsteady flow[END_REF] used leading edge actuation to recover lift on a 3D airfoil after separation of the boundary layer. Colonius and Williams [START_REF] Colonius | Control of vortex shedding on two-and threedimensional aerofoils[END_REF] puts forward that the dynamics response of lift shows some universal features, whatever the actuation or the type of profile of interest [START_REF] Woo | Transitory control of a pitching airfoil using pulse combustion actuation[END_REF][START_REF] Williams | Closed-loop control of a wing in an unsteady flow[END_REF]. They observed a response time of about a dozen of chord convective time t * = C/U ∞ , where C is the chord and U ∞ is the freestream velocity. The same order of magnitude is also observed by Darabi and Wygnanski [START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 1. the forced reattachment process[END_REF][START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 2. the separation process[END_REF] in a separated flow over a flap, by Mathis et al. [START_REF] Mathis | Experimental study of transient forced turbulent separation and reattachment on a bevelled trailing edge[END_REF] for a separation on a bevelled trailing edge and by Shaqarin et al. [START_REF] Shaqarin | Open and closed-loop experiments to identify the separated flow dynamics of a thick turbulent boundary layer[END_REF] in a thick separated boundary layer over a flap. Moreover, it is generally observed that a pulsed fluidic actuation has increased performances compared to a continuous fluidic actuation [START_REF] Jones | An active flow circulation controlled flap concept for general aviation aircraft applications[END_REF][START_REF] Shaqarin | Open and closed-loop experiments to identify the separated flow dynamics of a thick turbulent boundary layer[END_REF]. Initial vorticity when the jet is blown might explain the extra suction obtained by pulsed jets. Similar to the observation of initial lift gain of controlled stalled airfoils made by Colonius and Williams [START_REF] Colonius | Control of vortex shedding on two-and threedimensional aerofoils[END_REF]. It seems nevertheless that the research effort in closed-loop lift control has mainly been focused on post-stall lift recovery.

Recently, Troshin and Seifert [START_REF] Troshin | Performance recovery of a thick turbulent airfoil using a distributed closed-loop flow control system[END_REF] performed closed-loop lift control on an airfoil with distributed fluidic actuation on the suction side of an unstalled thick wing, aiming at the recovery of lift loss due to surface degradation. They observed the lift to respond as a second order system with a typical response time of 8.1 convective chord time (i.e. in the same order of magnitude of the stalled cases). They also showed that the system dynamics may vary with both the Reynolds number and the angle of attack (AOA). The higher is the AOA, the higher is the gain, with a possibly faster response when AOA or Reynolds number increases (it is not possible to distinguish AOA and Reynolds number effects from their results).

Therefore, the aim of this paper is to provide a steady-state as well a dynamical information of lift control of a wind turbine blade profile whose trailing

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
edge is modified to allow circulation control. The control objective is to compensate atmospheric perturbations from manipulation of the lift, and closed-loop control is considered to dynamically maintain this objective. The experimental apparatus is described in the section II, then the average natural and actuated flow fields are described using PIV results in section III, providing information on the actuation functioning. In section IV we provide steady state information on lift control. Section V is dedicated to the dynamics of lift control in both open-and closed-loop configurations. Concluding remarks are given in section VI.

Experimental set-up

Blade profile

The trailing edge of a N ACA-65 4 421 airfoil profile has been rounded (radius of curvature 2% of the chord) and the camber of the airfoil has been adjusted to compensate aerodynamic degradation due to the round trailing edge. The profile is presented in figure 1 (left). The standard NACA-65(4)421 profile is a thick profile with two C L drops corresponding to a first boundary layer separation at the trailing edge of the profile for AOA 8˚, called latter α sep , and a second strong flow separation at the leading edge for AOA 20˚ [START_REF] Devinant | Experimental study of wind-turbine airfoil aerodynamics in high turbulence[END_REF].

From α sep to 20˚the separation point moves gradually from the trailing edge to the leading edge, corresponding to a gradual variation of the loads. As can be seen in figure 1, this is also the case for the modified NACA-65(4)421 of the present work as shown by Balériola [START_REF] Balériola | Etude expérimentale dactionneurs plasma surfaciques pour optimiser les performances aérodynamiques dun profil éolien[END_REF] and Aubrun and Leroy [START_REF] Aubrun | Aerodynamic loads fluctuations alleviation by circulation control on wind turbine blades[END_REF]. The relation between the two C L drops and the flow separation state of the current profile has been investigated by Balériola [START_REF] Balériola | Etude expérimentale dactionneurs plasma surfaciques pour optimiser les performances aérodynamiques dun profil éolien[END_REF] using PIV measurements. We focus our wing.

The discrete microjets are distributed along the wing span from the stack of thin profiles (1 mm spanwise) and thick blind profiles (20mm spanwise). The thin profiles possess a round hole connected to a small air circuit (of 1 mm width) open to the surface of the profile. This is illustrated in figure 2 (Left).

The thick profile possess the same round hole with no connection to the surface.

Once stacked, the holes form a plenum chamber in the span direction of the profile, which can be connected to pressurized air to generate discrete squared fluidic microjets. The same principle is reversely used to measure the pressure distribution around the chord. The measured pressure is hence integrated along the span of the profile, providing results close to what can be obtained in a 2D pressure field.

Wind tunnel

The section is below 0.4% and the operating speed of the wind tunnel is set to

U ∞ = 10m/

Unsteady angle of attack

A step in Reynolds number would be pefectly suited to model a gust effect.

However, it is not trivial to perform a wind step with a standard wind tunnel.

At the blade scale, a simple steady shear inflow in front of a wind turbine rotor disk induces a cycling variation of the blade incidence [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. Therefore, some aerodynamic facilities have modeled mean wind shear inflows effects using an active rotation of the wing [START_REF] Melius | Dynamic stall of an experimental wind turbine blade[END_REF]. Also, additional devices can be used at the inlet of the test section to explore different slow harmonic perturbations [START_REF] Tang | Experiments and analysis for a gust generator in a wind tunnel[END_REF]. These studies do not however take into account the additional impact of small turbulence scales. Impact of free-stream turbulence has been investigated for the present NACA profile by Devinant et al. [START_REF] Devinant | Experimental study of wind-turbine airfoil aerodynamics in high turbulence[END_REF]. They show that high turbulence levels have a major effect on aerodynamic properties when the separation point moves from the trailing edge to the leading edge (AOA > 8˚). At the rotor scale, active grids are nowadays used to reproduce spectral features similar to those found in the atmospheric boundary layer [START_REF] Reinke | Flow modulation by an active grid[END_REF]. These perturbations can adequately model stationary mean and turbulent flow properties of atmospheric

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT turbulent flows.
A the time of writing this article, a perturbation system is under development in our wind tunnel facility to take into account the impact of more realistic gust conditions on the control efficiency. However, a first simplified approach is considered in the present work. The wing is mounted on a specific bench that automatically controls the AoA of the wing to simulate a constant shear inflow in front of the wind turbine rotor. We focus on the linear part of the lift curve (AOA < 8˚) with only modest explorations towards AoA for which the flow separates at the trailing edge (i.e. AOA 10˚but

AoA << 20˚).

The wing can oscillate at a chosen frequency and amplitude of the incidence, around a specified AOA. The maximum allowed frequency is 1 Hz (or St max = f C U = 0.03) and a maximum amplitude of 6˚was used in the present study.

Air delivery system

The compressed air circuit delivering the air to the actuator is composed of a pressure regulator, a plenum tank and the actuation is monitored through the help of electric valves. The pressure regulator is used to expand the incomming compressed air to a given pressure, monitoring the actuation intensity by setting the microjet stagnation pressure. The expanded air is then stocked in an intermediate plenum tank working as a decoupling capacitor, ensuring the stability of the actuator stagnation pressure while the flow is rapidly set when opening the valves. The momentum coefficient C µ is used as a measure of the actuation strength. It is defined as :

C µ = Q µj V µj 1/2ρ ∞ U 2 ∞ S ref , (1)
where, Q µj is the injection mass flow rate, V µj is the injected flow velocity,

S ref = C • L
is the reference area of the wing, and the subscript ∞ refers to the external flow conditions. It is straightforward to show that:

C µ = 2 • P µj M 2 µj P ∞ M 2 ∞ • S µj S ref , (2)

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

where P µj is the static pressure of the microjet at the exit section, M µj is the exit Mach number of the microjet and S µj the exit section area of the actuators.

The present experimental set-up allows to vary C µ from 0 to 0.14.

As long as the injection stays subsonic, P µj is equal to the pressure at the location of the actuator. Hence, at constant incoming flow conditions, fixed AOA and fixed microjet exit section, C µ solely depends on the ratio of microjet to external Mach numbers M µj /M ∞ . This ratio can only be changed by increasing the stagnation pressure of the micro-jet. The relationship between C µ and the stagnation pressure of the micro-jets is theoretically non-linear, but in the limited range of M µj that could be set in the present experiment C µ can be well approximated by a linear function of the microjet stagnation pressure.

Note that we used the compressible form of the C µ because the micro-jet Mach number can easily exceed M µj = 0.3 (the respective stagnation to static pressure ratio is 1.064). Furthermore, since M µj is directly obtained from the isentropic relations, there is no need to know the temperature of the microjet flow to compute the C µ , which makes this formulation easier to use.

Measurement apparatus 2.4.1. Balance

First, the profile was placed on a 6-axis balance system to measure the time averaged lift and drag. The balance was carefully calibrated so that lift and drag uncertainties were estimated to be less than 2% in the free stream velocity range used in this study. with this configuration. Note that this sampling frequency is fast enough with respect to the characteristic time scales of interest in this paper. At each rising edge of the trigger signal, the pressure scanner performed an acquisition of all the channels and sent the digital values over the network to a dedicated computer. The number of received samples were compared to the number of rising edges to control that no packet losses were encountered.

The closed loop controler we used was integrated in a PC running a Linux OS with a RT PREEMPT patch for hard real time performances. This computer was equipped with an ADC PCI board for the input-output purposes. The expected input of this controler is analog voltage values from pressure transducers.

Since the pressure scanners embed their own ADC, it was not possbile to use them in our feedback-loop. Therefore a set of 16 FirstSensor analog pressure transducers was used, allowing to measure the pressure using a standard acquisition board on a PC. These sensors have a response frequency of 2kHz, sligthly faster than the pressure scanners.

PIV

The flow field around the actuator was probed using Particle Image Velocimetry (PIV). A 2D-2C PIV system was used for cross-sections of the profile 6 µm slightly smaller than the pixel size [START_REF] Adrian | Particle-imaging techniques for experimental fluid mechanics[END_REF]. Nevertheless, thanks to the sufficient amount of light provided by the 200 mJ laser, the actual particle images were bigger and no sign of peak-locking could be observed in the data. The flow was seeded with olive oil particles generated through Lazkin nozzles. Their averaged diameter is less than 1 µm giving a response time of 2 µs which is sufficient for the time scales of interest in this paper. We set the inter-frame time accordingly to the maximum displacement and/or out of plane motion to maximize both the measurement accuracy and the cross-correlation value. The inter-frame time ranged from 5 to 20µs.

For the Stereo-PIV the lenses were mounted to the cameras with Scheimpflug adaptors. It was possible to set both camera at a 45 degree viewing angle, ensuring that measurement errors are well balanced between all the velocity components [START_REF] Willert | Stereoscopic digital particle image velocimetry for application in wind tunnel flow s[END_REF][START_REF] Prasad | Stereoscopic particle image velocimetry[END_REF].

The PIV image processing was performed using an iterative algorithm with deforming interrogation windows [START_REF] Huang | On errors of digital particle image velocimetry[END_REF][START_REF] Scarano | Iterative image deformation methods in piv[END_REF] and seeding adaptivity [START_REF] Wieneke | Adaptive piv with variable interrogation window size and shape[END_REF]. The initial window size was set at 32x32 pixels with an overlap of 50%. Outliers were rejected using a correlation peak-ratio of 1.3 together with a Universal Outlier Detection spatial filter [START_REF] Westerweel | Universal outlier detection for piv data[END_REF]. The processing was carried out in both Davis expands for about 5% of the chord length downstream, which represent about one trailing edge diameter and is a typical bluff body recirculation bubble size.

When the control is on, the presence of the micro-jet flow is evident from the yellow trace that is left in the streamwise velocity iso-contours. Contrary to what was originally expected, the micro jet shows a limited turning angle 265 around the trailing edge Coanda surface: it separates rather early from the wall. This is due to a combination of a too high microjet height relative to the Coanda surface radius, and three dimensional effects from the finite span of the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

discrete actuation slots [START_REF] Kadosch | Déviation des jets par adhérence à une paroi convexe[END_REF][START_REF] Englar | Experimental investigation of the high velocity coanda wall jet applied to bluff trailing edge circulation control airfoils[END_REF]. Nevertheless, we can see from figure 4 that the actuation diminishes the overall size of the recirculation bubble and displaces the wake towards the high pressure side of the wing, which is highlighted from the change of sign of the vertical velocity component in the wake (see figure 4 (Right)). Therefore, even if the actuation do not seem to perform as initially expected and that the circulation control can probably be improved using a better suited microjet/Coanda-surface pair geometry, the actuation solution proposed here has a strong effect on the flow and it will be used later to analyze the dynamic response of the flow over the blade for pre-stall configurations (see section 5).

Cross-flow velocity field

As presented in the section 2 the fluidic actuation is three dimensional in our study, hence we need to verify how it performs along the span and seek for possible three dimensional effects. To summarize, we have shown that the proposed circulation control system affects the averaged flow field by reducing and by imposing a three-dimensional effect on the recirculation area behind the trailing edge. It has also been shown 305 that these effect are localized in the neighboring area of the microjet main axis.

In the following, we want to provide a quantification of the actuation effect on From this figure, we can see that the actuation has numerous effects on the lift curves. First of all, as expected from previous circulation control attempt, the lift is increased for all angles of attack [START_REF] Seifert | Oscillatory blowing: a tool to delay boundary-layer separation[END_REF][START_REF] Jones | An active flow circulation controlled flap concept for general aviation aircraft applications[END_REF]. This is an important characteristic if one aims at controlling the lift perturbations created by external wind gusts for example. It is also clear from the C L curves that the actuation plays a significant role in progressively delaying the threshold angle at which the onset of boundary layer separation, α sep , occurs at the trailing edge [START_REF] Balériola | Etude expérimentale dactionneurs plasma surfaciques pour optimiser les performances aérodynamiques dun profil éolien[END_REF]. It should be noted however that it was not possible with the available blade design to blow jets that decrease lift. A new design is under progress for this purpose and only Thrust compensated drag (left)

alleviation of gust that increase wind speed are considered here.

Another observation that can be made from figure 6 is that the static gain is not only dependent on the actuation strength C µ but also on the flow state (i.e. the angle of attack). As can be seen, much more lift is gained after the threshold angle α sep than at attached flow regimes (α ≤ α sep ). This behavior was already observed for various CC devices [START_REF] Greenblatt | The control of flow separation by periodic excitation[END_REF], therefore showing that the currrent CC strategy works with the same principles as what can be found in the litterature.

Note however that the maximum gain seems to be located around α = 10˚and reduces for higher AOA.

In figure 7 left we present a comparison of the C D curves for different C µ .

At a first glance, one could think that the control decreases the drag of the profile. However, this is a simple microjet thrust effect. Indeed, each individual microjet acts as a small nozzle that creates thrust in the opposite direction of drag, that can be measured with zero free stream velocity. By removing this thrust effect on the drag measurements, drag curves overlap for all controlled cases (see figure 7

Sensitivity of actuation

We plot in figure 8 the lift coefficient gain against C µ for various AOA in the attached regime (i.e. α ≤ α sep ). It is clear that the control gives a linear change in lift for all presented AOA. As stated earlier, the relationship between the micro-jet stagnation pressure and C µ can be well approximated by a linear relation. Since the relation between C µ and the micro-jet generating pressure is almost linear (see dedicated previous section), it means that the lift increment due to the actuation is directly proportional to the pressure set in the micro-jet stagnation chamber. This is an interesting consideration for a future feedback control strategy because pressure is a quantity easily measured.

We present on the same figures linear regressions of the measured lift increments. We provide the estimated slope ∂∆C L ∂Cµ in the legends of figure 8 and we see that the actuation efficiency is greater for higher angle of attack (from left to right). Note that this trend was also obsered by Troshin and Seifert [START_REF] Troshin | Performance recovery of a thick turbulent airfoil using a distributed closed-loop flow control system[END_REF] in a different control configuration. In our case, the natural position of the stagnation point at the trailing edge moves upstream with increasing AOA. This provides more room for the control to displace it downstream, allowing for more circulation changes at higher AOA.

Discussion on the performance of the current actuation

The comparison of performances with literature is a difficult exercise because of the variety of the different techniques that were used to manipulate

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

lift. In a review, Greenblatt and Wygnanski [START_REF] Greenblatt | The control of flow separation by periodic excitation[END_REF] proposed two main parameters to evaluate the performances of various actuation devices proposed in the literature. The first one ∆C L,max corresponds to the gain in maximum lift, and the second one is the maximum post-stall lift gain ∆C L . We will only use the first one here, since we do not deal with stalled cases in our current study. For example, Seifert et al [START_REF] Seifert | Oscillatory blowing: a tool to delay boundary-layer separation[END_REF] reported a 21% increase of ∆C L,max for a C µ = 0.001 using a 2D continuous injection on a flapped NACA 0015 airfoil. From figure 6 we can see that the current actuation device offers a ∆C L,max of 32% but for a C µ = 0.139. Additionnaly, Cagle and Jones [START_REF] Cagle | A wind tunnel model to explore unsteady circulation control for general aviation applications[END_REF] and Jones et al. [START_REF] Jones | An active flow circulation controlled flap concept for general aviation aircraft applications[END_REF] used the ratio of lift gain over injection coefficient ∆C L ∆Cµ to evaluate the actuation efficiency at a given AOA. They reported a ratio of 50 at α = 0˚, for a similar concept of lift control but using 2D actuation. Here, we measured a more modest but non negligible ∆C L ∆Cµ of 1.29 at the same AOA. It is believed that these discrepancies are due to the three-dimensional actuation performed in our study. Let's assume that the circulation control is localized in span along the 2D sections of the foil where the micro-jets exists, which is not a strong assumption given the observations made in the previous section.

We can therefore recompute the lift gain using the sum of the individual orifice spans, l a , instead of the wing span, L, as follows: ∆C L,max2D = ∆C L,max • L/l a . This leads to a gain of 840% of ∆C L,max2D at C µ = 0.139 corresponding to a non negligible 6% increase of ∆C L,max2D for C µ = 0.001. This is closer to the reported gains of Seifert et al. [START_REF] Seifert | Oscillatory blowing: a tool to delay boundary-layer separation[END_REF]. Similarly, we then have a 2D lift to momentum coeficient ratio of ∆C L,2D ∆Cµ = 26, which is of the same order of magnitude than the results of Cagle and Jones [START_REF] Cagle | A wind tunnel model to explore unsteady circulation control for general aviation applications[END_REF].

This shows that the chosen actuation method performs in a comparable manner as the different techniques reported in the literature and confirms that most of the lift gain is actually obtained through local modifications of the circulation, as we inferred from the observations of the PIV results in the previous section. x/C = 0.09

M A N U S C R I P T A C C E P T E D

P(t)-P i [Pa] t*=tU/C [s] measured -5
.67 e -t/3.46e+00 x/C = 0.89

P(t)-P i [Pa] t*=tU/C [s] measured -8
.88 e -t/5.54e-01 As can be seen, the pressure response for detached flow cases no longer resembles to a response of a first order ODE system. At these AOA, it is possible to recognize an under-dumped second order ODE system response, with a small overshoot.

In order to have a better overview of the multiple outputs of the system, we x/C = 0.09

P(t)-P i [Pa] t*=tU/C [s] measured -14
.84 e -t/2.08e+00 x/C = 0.89

P(t)-P i [Pa] t*=tU/C [s]
measured -14.18 e -t/1.33e+00 fitted the following first order type response for each pressure response:

∆P (x, t) = P (x, t) -P i = K u(t) 1 -e -t/τ , (3)
where P is the pressure, t the time, K = K(x) the steady state gain, τ = τ (x) the time constant of the system (i.e. time to reach 63% of the steady state), P i = P i (x) the initial value of the pressure over the upper surface of the blade and u(t) is the control input that is equal to unity when actuation is ON and to zero otherwise. Typical fitted results are given in figures 9 and 10. Even though, this response type is not adapted to the pressure response at high angles of attack, it can be seen in figure 10 that the error between the fitted function and the measured response is relatively small such that the gains and response times obtained are fairly representative of the real system. Results can be separated in two areas. The first one near the actuator were a significant peak, proportional to the control input C µ , is observed and should induce an increase of the pitching moment as usually observed in CC devices [START_REF] Englar | 2-d circulation control airfoil benchmark experiments intended for cfd code validation[END_REF]. The second area is plotted in figure 11 using the different dimensionless steady state gains found around the chord of the profile for different control inputs (i.e. C µ ) at α = 0.

First of all, contrary to what is reported by many authors [START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 1. the forced reattachment process[END_REF][START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 2. the separation process[END_REF][START_REF] Mathis | Experimental study of transient forced turbulent separation and reattachment on a bevelled trailing edge[END_REF], we didn't find significant differences between dynamic responses to an actuation step ON -K/(0.5 ρ (open symbols of figure 11) and an actuation step OFF (closed symbols of figure 11). The main difference between these studies is the initial state of the flow before control. In the results of figure 11 the flow is attached all the time while in the other studies the flow is either initially separated [START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 1. the forced reattachment process[END_REF][START_REF] Darabi | Active management of naturally separated flow over a solid surface. part 2. the separation process[END_REF], or the final targeted state is the flow separation [START_REF] Mathis | Experimental study of transient forced turbulent separation and reattachment on a bevelled trailing edge[END_REF].

U 2 ) x/C C µ =0.051 C µ =0.112 C µ =0.139 C µ =0.051 C µ =0.112 C µ =0.139
The second important observation is a clear effect of the actuation intensity on the steady state gains: the stronger the actuation, the higher the steady state gain. This is in agreement with the previous results of the balance measurements (see figure 8). At last, results show that the respective contributions of both areas (i.e. nearby the actuation and around the chord) to the integrated gain (lift coefficient) are of the same order of magnitude. Therefore, the actuation used in this study, although very local, modifies the overall pressure field around the wing, similarly as can be found in other circulation control configurations [START_REF] Jones | An active flow circulation controlled flap concept for general aviation aircraft applications[END_REF][START_REF] Cagle | A wind tunnel model to explore unsteady circulation control for general aviation applications[END_REF]. time needed to reach 63% of the steady state. The 99% steady state response time, τ 99 , is three times τ . τ 99 is found to range from 1.5 to 13.5 convective times. Note that previous studies did not investigate the response time around the chord but rather close to the actuation, however, for locations close the leading edge (or away from the actuation) the response time of pressure is of the same order of magnitude of the one reported in the literature in separated flow over a wing [START_REF] Colonius | Control of vortex shedding on two-and threedimensional aerofoils[END_REF][START_REF] Williams | Lift response of a stalled wing to pulsatile disturbances[END_REF] and in detached boundary layer flow [START_REF] Mathis | Experimental study of transient forced turbulent separation and reattachment on a bevelled trailing edge[END_REF][START_REF] Shaqarin | Open and closed-loop experiments to identify the separated flow dynamics of a thick turbulent boundary layer[END_REF].

Results of the present study show a clear trend of the response time around the chord. For all input values, the time constant is increasing with the distance from the actuation. If the static gains were homogeneously distributed around the profile, no steady states change in pitching moment is expected. However, the fact that the pressure takes longer to establish far from the actuation indicates that transient pitch torque is created during the establishment of the modified circulation. This is rather important behavior for final application, that might also be found in other CC devices.

We present in figure 13 the chordwise distribution of K and τ measured with C µ = 0.139 for varying AOA. As expected for the lift curve presented in figure 6, the gain depends highly on the AOA on the suction side, while almost no effect is observed on the pressure side. The highest gain is found for α = 10 in good agreement with the gain of C L found in figure 6. For higher AOA, the gain deacreases, also in good agreement with figure 6.

Regarding the response time of pressure, no effect of the angle of attack on the response time is perceived for attached flow cases. As previously observed, longer response times are measured for locations farther from the actuator. This observation is no longer valid for detach regimes where the response time seems to be more evenly distributed along the chord. The α = 10˚results show that for x/C < 0.6 the response times are comparable to smaller AOA whereas for x/C > 0.6 the response times are increased compared to the attached flow cases.

For α = 14, the response times are clearly slowed near the trailing ege (i.e. for x/C > 0.5) and accelerated near the leading edge (i.e. for x/C > 0.5) inducing a relativelly constant time response over the chord. This is an indication that the separated flow possesses a different dynamics than the attached one.

The results of this section show that the dynamic response of the flow over an airfoil is strongly dependent on the presence of flow detachment on the suction side of the wing. Even though no direct measurement of the pitching moment was performed, the results of this section show that a transient pitching mo-

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

ment is expected from the local response time of the pressure that depends on the distance to the actuation. Hence, the fact that the local response time of pressure depends on the distance to the actuation is an interesting result that should be beard in mind.

Nevertheless, it seems that the pressure response to the actuation can be reasonably approximated by typical first order system types, demonstrating that the dynamics of pressure establishment remains relatively simple.

Phase locked C L response

From the former phase locked pressure step responses, it is possible to build the corresponding lift step responses by simply integrating the pressure field on the wing surface :

C L (t) = 1 S ref S C p (x, t)dx (4)
where

C p = P -P∞ 0.5ρU 2 ∞
the pressure coefficient with x the projection of the profile skin coordinates in the streamwise direction. The lift response ∆C L (t) can then be directly obtained by spatially integrating the measured pressure response at each time step:

∆C L (t) = 2 ρU 2 ∞ S ref S ∆P (x, t) dx, (5)
where ∆P (x, t) is defined in equation 3.

The measured lift step responses are plotted in figure 14 for the three different momentum coefficients at α = 0. The lift gain is normalized by the steady state gain and all the lift responses collapse onto a single curve corresponding to a typical step response from an ODE system. A least-square approximation of equation C L (t) = 1 -K L e -t/τ L is performed and plotted against experimental measurement points. The response time to reach 63% of the steady state is found of the order of one convective time τ = U∞ C while the response time to reach 99% of the steady state is more than 3 convective times.

Theoretical lift response to gusts for attached flow cases is usually given by the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Küssner's function [START_REF] Küssner | Zusammenfassender bericht über den instationären auftrieb von flügeln[END_REF] obtained using potential flow theory. The derived lift response can be approximated by a three exponential term function [START_REF] Jones | The unsteady lift of a wing of finite aspect ratio[END_REF]:

C Lg (s) = 1.0 -0.236e 0.058s -0.513e 0.364s -0.171e 2.42s , (6)
where s = t C 2U∞ is the dimensionless time based on the freestream velocity and the half chord. The former equation is also reported on figure 14.

As can be seen, the measured lift response is faster than what is theoretically predicted. There might be several reasons why the theory fails in recovering For detached cases, the response curves are different and show an overshoot of lift before it reaches the steady state value, which is a typical behavior of under-damped second order system. The response time of lift at these angles of attack is therefore slightly longer than for attached cases. As a matter of comparison, we derived the step responses from the indicial lift response measured by Colonius and Williams [START_REF] Colonius | Control of vortex shedding on two-and threedimensional aerofoils[END_REF] on a semi-circular wing with leading edge flow separation. This curve is reproduced on figure 15 and also exhibits an overshoot but with a slower response than for our cases for which only trailing edge flow separation occurs.

These result may suggest that the lift response strongly depends on whether the flow is separated or attached but also on how/where it separates from the wall. Separation close to the trailing edge seems to provide faster responses than separation at the leading edge of the wing.

Comparison with turbulent scales

The response time scale observed in the previous paragraph has to be compared with the times scales of perturbations that can be found in the atmospheric boundary layer. Typical spectrum of atmospheric boundary layer velocity fluctuations are given in figure 16, based on the formula given by Kaimal et al [START_REF] Kaimal | Spectral characteristics of surface-layer turbulence[END_REF]. This spectrum is plotted against a reduced frequency representative of the characteristic scales of an existing wind turbine blade with chord length of 1m and velocity of 60m/s. As can be seen, most of the velocity perturbations are low frequency with respect to the blade time scales: the peak frequency of all the perturbations are found below St < 0.1 for all the 3 velocity components. Since the open-loop control system has a characteristic frequency of St = 1, it is clear that the proposed system is rapid enough to account for those scales. This is illustrated in figure 16, where the shaded area corresponds to frequencies above the cut-off frequency of the open-loop system, and the transfer function of the lift system is presented in dark-grey. It should be mentioned that the response time of closed-loop systems can be made faster than the open-loop one, hence ensuring that it is physically possible the account for external perturbations at the blade scale.

Closed-loop control

Control Modeling

From the results of the former section a numerical model of the system can be built to evaluate the performances of a feedback controller. Since the system can be represented by a first order ODE system, a PID controller should be able to handle the system [START_REF] Shaqarin | Open and closed-loop experiments to identify the separated flow dynamics of a thick turbulent boundary layer[END_REF]. C Lp = 0.05 cos(w p t * ) with w p the perturbation frequency chosen according to figure [START_REF] Shires | Application of circulation controlled blades for vertical axis wind turbines[END_REF]. As can be seen, even with a binarized actuation the actuator is still able to minimize the effect of the perturbations so that the output lift stays close to the setpoint.

It is clear form figure 18 that the control could not entirely remove the lift perturbation. The remaining lift oscillations are due to the non-proportional behavior of the valves, indeed if the valves are modeled as ideal ones (i.e. have an infinitesimal reaction time), the incoming perturbations could be entirely removed as can be seen in figure 19.

Another important point that arises from the fact that our actuators are binary, is that parameters of the controller (i.e. PID) have very limited impact on the efficiency of the control. This is obvious considering that the actuation signal computed by the PID is actually binarized a posteriori by the valves. This Since the actuation can only increase the lift, the perturbation consists in a sinusoidal variation of α defined as follows:

α(t) = α 0 -∆α • sin(2πf α t), (7)
where α 0 is the AOA at which the lift setpoint is measured, ∆α is the amplitude of the perturbation and f α the oscillation frequency. The oscillation frequency is chosen so that St α = fαC U = 0.003. It is at the same time low enough to avoid dynamic stall effect and/or hysteresis effect that can be expected to appear around St = fαC U = 0.1 [START_REF] Choudhry | An insight into the dynamic stall lift characteristics[END_REF][START_REF] Mccroskey | An experimental study of dynamic stall on advanced airfoil sections[END_REF][START_REF] Mcalister | An experimental study of dynamic stall on advanced airfoil section[END_REF], but high enough to be representative of scaled atmospheric turbulent perturbations.

The feedback loop used in this experiment is shown in figure 20 (left). The error between measured lift and the setpoint is computed and fed to a proportional controller. The output of the controller is then binarized to a TTL electric signal that can be used to actuate the electric valves. Then a new estimate of lift is computed by integrating the pressure on the wing and compared to the setpoint, hence closing the loop. The results of the closed-loop control in the attached regime are presented in figure 21 for α sp = 4˚and ∆α = 2˚. Before the control is turned on, the lift oscillates in phase with the perturbation. It is clear that, as soon as the control is ON, the lift oscillations are reduced by 83%, proving that feedback control of lift using fluidic actuation is an efficient method. The results of the

PI

L 0 C µ - + L C µ * ΔL L(t) U ∞ α(t)

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

closed-loop control in the separated regime is presented in figure 22. The main difference with the previous results is that there is a phase opposition between the AOA perturbation and the lift output before control. This is expected from the change of slope in the C L vs. α curve for AOA greater than 10˚. However, this does not modify the control efficiency which is found of the same order of magnitude. These results not only show the feasibility of feedback lift control but also show that this simple control can be applied even when the regime of the flow over the blade is in the non-linear part of the lift coefficient.

Conclusions and Perspectives

5 1 M

 1 loads on the turbine that can lead to severe damages. Controlling the lift of the individual wind turbine blades may thus be considered to soften these undesired

Figure 1 :

 1 Figure 1: (Right) Arifoil profile. The original N ACA -65 4 421 airfoil is in red dash line and the control profile N ACA -65 4 421 -CC, used in this study, is in solid black line. (Left) modified profile aerodynamic performances from Aubrun and Leroy [33].

Figure 2 :

 2 Figure 2: Left: diagram of a thin (non-blind) profile. Right:photograph of the thick (blind) and thin (non-blind) profiles.

  s in the present study. The 2D blade is mounted between two vertical flat planes in order to achieve a 2D flow configuration. The wing span was equal to the distance between these two plates and measured L = 1.1m. At this free stream velocity, the Reynolds number of the flow based on the chord length (C = 0.3 m) is slightly above 2 × 10 5 .

2. 4 . 2 .

 42 Wall pressure sensing Two sets of pressure transducers were used during the experiments. The first one was used during the open loop characterization of the lift dynamics. The pressure taps were then connected to two Chell MicroDAQ pressure scanners to synchronously obtain the pressure field around the profile. The scanners were externally triggered at 500Hz, the maximum possible sampling frequency M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

(see figure 3 (Figure 3 :

 33 figure 3 (Right)). In both cases the cameras and the laser were placed on a traverse system allowing to change the measurement position without modify the calibration.

8 . 2 from 3 .Figure 4 :

 8234 Figure 4: Isocontours of the averaged streamwise (left column) and vertical (right column) velocity components for various injection coefficients: Cµ = 0 top row, Cµ = 0.051 middle row, Cµ = 0.11 bottom row. Results obtained from 2D-2C PIV in the center plane of the wing z = 0 (i.e. including an actuation micro-jet) and A0A = 0. The black line represents the U = 0 iso-contour, delimiting the recirculation bubble.

Figure 5 :

 5 Figure 5: Cross-flow velocity fields taken at x/C = 1.03 and AOA = 0 for the natural (left column) and for the actuated Cµ = 0.11 (right column) cases. The iso-contour maps represent the streamwise (top row), vertical (middle row) and spanwise (bottom row) velocity component.

Figure 6 : 4 .

 64 Figure 6: Comparison of actuated versus non-actuated lift curves for the current wing. (left) 10m/s

Figure 7 :

 7 Figure 7: Comparison of actuated versus non-actuated drag curves for the current wing (left).

Figure 8 :

 8 Figure 8: Lift coefficient as function of injection coefficient for α = 0 (left), α = 5 (center), α = 8 (right).

Figure 9 :

 9 Figure 9: Pressure step response to actuation on the upper side of the wing at x/C = 0.09 (left) x/C = 0.89 (right). The results are presented for Cµ = 0.139 and α = 0.

Figure 10 :

 10 Figure 10: Pressure step response to actuation on the upper side of the wing at x/C = 0.09 (left) x/C = 0.89 (right). The results are presented for Cµ = 0.139 and α = 14.

Figure 11 :

 11 Figure 11: Evolution of the gains K with the chord for varying actuation intensities, the vertical axis was chose to cut-off the very high gains obtained in the neighborhood of the actuator, open and closed symbols are for opening and closing actuation respectively

Figure 12 showsFigure 12 :

 1212 Figure 12 shows the dimensionless time constants τ found around the chord of the profile for different control inputs. This response time corresponds to the

Figure 13 :

 13 Figure 13: Left: chordwise evolution of the gain K for different angles of attack. Right: evolution of the response time τ with the chord for different angles of attack. On both figure attached regimes are presented with black symbols, whereas blue symbols represent detached flow regimes.

  the measured response time. First of all the theoretical model uses inviscid flow and thus cannot model the fact that the CC actuation modifies the flow separation on the trailing edge. Moreover, the theory and the experiment do not use identical type of perturbation: in the theoretical model the gust is modeled as a local camber modification of the wing whereas in the experiment the perturbation is introduced at the trailing edge. It is therefore difficult to conclude on the observed differences but it must be noticed that blade profiles lift might have a faster response time than what is generally considered for fluid/structure interaction.

Figure 14 :

 14 Figure 14: Step responses of the lift coefficient C L for different actuation intensity obtained by integrating all the measured pressure responses at α = 0.

Figure 15 :

 15 Figure 15: Step responses of the lift coefficient C L obtained by integrating all the measured pressure responses for various AOA at Cµ = 0.14 compared with theoretical functions as well as other experimental results. For sake of clarity only one point over 10 is plotted.

Figure 16 :Figure 17 :

 1617 Figure 16: Spectrum velocity perturbations in an atmospheric boundary layer taken at an altitude of z = 80m and average velocity of 8m/s. The frequency is reduced by a wind turbine blade length and velocity scales. The grey curve represents the efficiency of the openloop control system: beyond St = 1 the control system becomes inefficient.

  wing at a given angle of attack, this estimation of the lift information (using the suction side only) was found enough for an efficient and simple proportional controller implementation.

Figure 20 :

 20 Figure 20: Diagram of the feedback-loop used on the experiment (left), sketch of the external lift perturbation applied to the wing by varying the angle of attack (right).

Figure 21 :

 21 Figure 21: Time traces of normalized lift and angle of attack. The lift controller is set at t = 40s. Setpoint measured at α = 4˚: the flow is attached.

Figure 22 :

 22 Figure 22: Time traces of normalized lift and angle of attack. The lift controller is set at t = 40s. Setpoint measured at α = 12˚: the flow is separated.

Acknowledgments

The authors want to thank the French National Research Agency for supporting this study through the SMARTEOLE project (ANR-14-CE05-0034). Also, the authors would like to thank the PRISME laboratory staff where the wind tunnel experiments presented here were performed.

means that no fine tuning of the controller parameters is actually needed when electromagnetic valves are used.

Experimental results

Based on the encouraging results obtained using the model of the previous 600 section, an experiment has been designed to test the closed-loop control. The wing is mounted on a motorized bench allowing to vary the angle of attack, which acts as a perturbation in lift as sketched in the figure 20 (right). The profile is equipped with 15 pressure sensors on the suction side of the wing which serve to estimate the lift. Note that no sensor was placed on the pressure side of


domain: spi phys

text:
Thomas Hainaut
email: thomas.hainaut@onera.fr.

Thomas Le Garrec
email: thomas.legarrec@onera.fr.

Cyril Polacsek
email: cyril.polacsek@onera.fr.

Daniel Mincu
email: daniel-ciprian.mincu@onera.fr.

Sébastien Deck
email: sebastien.deck@onera.fr.

Aerodynamic and Aeroacoustic Numerical Investigation of Turbofan Engines using Lattice Boltzmann Methods

come   L'archive ouverte pluridisciplinaire

I. Introduction

With the ever increasing air traffic and stricter noise regulations around airports, aircraft manufacturers require faster and more accurate numerical tools to allow for reliable aerodynamic and aeroacoustic investigations of complex systems. In recent years, Lattice Boltzmann Methods (LBM) proved efficient at simulating very large CFD problems on massively parallel clusters [START_REF] Khorrami | Towards full aircraft airframe noise prediction: lattice Boltzmann simulations[END_REF][START_REF] Sengissen | Simulations of LAGOON landing-gear noise using Lattice Boltzmann Solver[END_REF][START_REF] Romani | APPLICATION OF LATTICE-BOLTZMANN METHOD TO ROTORCRAFT AERODY-NAMICS AND AEROACOUSTICS[END_REF][START_REF] Barad | Lattice Boltzmann and Navier-Stokes Cartesian CFD Approaches for Airframe Noise Predictions[END_REF][START_REF] Shin | Aircraft Cabin Outflow Valve Tonal and Broadband Noise Prediction Using the Lattice Boltzmann Method[END_REF]. One advantage of LBM is its natural very low dissipative method, allowing for relatively small number of points per wavelength. However, one important drawback, which makes it strength on other aspects, is its lattice structure which is, in itself, not perfectly adapted to account for complex geometries and mesh turbulent boundary layers. To solves this first issue, immersed boundary conditions are commonly employed [START_REF] Khorrami | Towards full aircraft airframe noise prediction: lattice Boltzmann simulations[END_REF][START_REF] Sengissen | Simulations of LAGOON landing-gear noise using Lattice Boltzmann Solver[END_REF][START_REF] Romani | APPLICATION OF LATTICE-BOLTZMANN METHOD TO ROTORCRAFT AERODY-NAMICS AND AEROACOUSTICS[END_REF]. To address this second issue, a turbulence model is performed in this work. It relies on a Delayed Detached Eddy Simulation (DDES) [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] based on a Spalart-Allmaras turbulence model [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF][START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF]. This turbulence model has been historically introduced for RANS solvers with the objective of computing the Reynolds stresses, as the averaged equations are not closed. Based on the Boussinesq hypothesis which states that he momentum transfer caused by turbulent eddies can be modeled with an eddy viscosity ν t , the Spalart-Allmaras model is a one transport equation for a viscosity-like variable, directly linked to the turbulent eddy viscosity ν t [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF][START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF].

In this work, in the context of turbomachine aerodynamic and aeroacoustic numerical investigations, a turbulence modeling at solid boundaries using a DDES approach is performed on three benchmark cases. Firstly, to assess the capability of self-noise generation, a symmetric aerofoil at high Reynolds number and at angle of attack is studied. Secondly, for leading edge noise predictions, a rod-aerofoil interaction noise case, where the rod generates a turbulence wake, which impinges the leading edge of the aerofoil. Also, the complex physics associated with the turbulence generation in the wake of the rod as well as the aeroacoustics from both the back of the rod and the leading edge of the aerofoil makes it a comprehensive benchmark. Finally, a case combining self-noise and interaction noise on mobile surfaces is investigated using a turbofan in rotation interacting with a grid-generated incoming turbulence.

II. Model

A. Boltzmann equation

In a standard CFD solver where the Navier-Stokes equations are solved, the governing equations describe the evolution of macroscopic values (pressure, density, velocity, temperature, ...), whereas in the Boltzmann equation, a mesoscopic description (kinetic energy, momentum, probability of collision between particles) is used. The general governing equation writes

∂ f (x, t, c) ∂t + c i ∂ f (x, t, c) ∂ x i = ∂ f (x, t, c) ∂t coll , (1)
where f (x, t, c) is the distribution function which represents the density of particles at point x, time t, with a velocity c.

The collision term, denoted as (∂ f (x, t, c)/∂t) coll in Equation ( 1), represents the rate of variation which redistributes the energy between particles. From the H-theorem [START_REF] Boltzmann | Weirere Studien uber das warmegleich-gewich unter gasmolekulen[END_REF], one can see that the collision term describes the tendency of the distribution function f (x, t, c) to tend to an equilibrium distribution function f eq (x, t, c), characterized by the Maxwell-Boltzmann distribution

f eq (x, t, c) = ρ m 2πk B T D/2 exp -m|c -u| 2 2k B T , (2)
where k B is the Boltzmann constant, T the temperature, m is the particle mass and D the degree of freedom of the particle.

The Bhatnagar-Gross-Krook (BGK) [? ] collision term is a popular [START_REF] Khorrami | Towards full aircraft airframe noise prediction: lattice Boltzmann simulations[END_REF][START_REF] Sengissen | Simulations of LAGOON landing-gear noise using Lattice Boltzmann Solver[END_REF][START_REF] Romani | APPLICATION OF LATTICE-BOLTZMANN METHOD TO ROTORCRAFT AERODY-NAMICS AND AEROACOUSTICS[END_REF][START_REF] Barad | Lattice Boltzmann and Navier-Stokes Cartesian CFD Approaches for Airframe Noise Predictions[END_REF][START_REF] Shin | Aircraft Cabin Outflow Valve Tonal and Broadband Noise Prediction Using the Lattice Boltzmann Method[END_REF] simplified linear approximation of the collision term which states that f (x, t, c) tends to f eq (x, t, c) in a characteristic time τ, often referred to as relaxation time. It reads

∂ f (x, t, c) ∂t coll = - f (x, t, c) -f eq (x, t, c) τ . (3)
To improve the robustness and accuracy of the BGK collision term, a variant based on a two-relaxation-time collision term [START_REF] Ricot | Lattice Boltzmann method with selective viscosity filter[END_REF] is implemented in this work. Referred to as Dual Relaxation Time (DRT), this method over-relaxes the distribution functions f (x, t, c) towards equilibrium while preserving the conservation of mass and momentum, and to reconstruct the non-equilibrium part through a second-order regularization process [START_REF] Ricot | Lattice Boltzmann method with selective viscosity filter[END_REF][START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF]. From a Chapman-Enskog expansion, one can retrieve the mass density ρ and the fluid momentum ρu from ρ(x, t) = ∫

B. Lattice Boltzmann method

In the present work, the Boltzmann equation is dicretized in 3D in a lattice structure along 19 discrete velocities c α [START_REF] He | Discrete Boltzmann equation model for nonideal gases[END_REF], as represented in Figure 1. Along each discrete velocity, Equation (1) can be seen as a first-order transport equation with a constant convection velocity. With such discrete velocities, the equilibrium distribution functions f eq α (x, t) in each α directions can be obtained using Hermite polynomials [START_REF] He | Discrete Boltzmann equation model for nonideal gases[END_REF][START_REF] Grad | Note on N-dimensional hermite polynomials[END_REF][START_REF] Grad | On the kinetic theory of rarefied gases[END_REF]

f eq α (x, t) = ρω α 1 + c α • u rT + 1 2 c α • u rT 2 - |u| 2 2rT , (5)
with ω α are the weighted constants associated with the Hermite polynomials and r = k B /m is the gas constant. Nevertheless, Equation ( 5) is truncated to the second-order, which introduces a approximation error in o(M 3 ), where M is the Mach number, which limits such discretization to Mach number up to 0.4.

For each discrete velocity c α , the associated density of particles is denotes f α (x, t), therefore equations ( 4) can be written

ρ(x, t) = α f α (x, t) and ρ(x, t)u(x, t) = α c α f α (x, t). (6)

C. Spalart-Allmaras turbulence model

To correctly simulate the turbulent eddy viscosity ν t , a fine grid is required near solid boundaries, however, with meshes composed of regular hexahedron elements, it quickly becomes too computationally expensive. To overcome this difficulty, the eddy viscosity can be modeled using a turbulence model, reducing the requirements on the mesh size. In this work, a DDES model [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] based on the Spalart-Allmaras turbulence model [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF][START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] is implemented. This model is originally based on the Spalart-Allmaras [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF][START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] RANS model which solves one transport equation for a eddy viscosity like variable ν. For clarity, the transition terms are omitted and the standard model reads

D ν Dt = c b1 S ν production + 1 σ ∆ • ((ν + ν)∆ ν) + c b2 (∆ ν) 2 diffusion -c ω1 f ω ν d ω 2 destruction . (7)
The turbulent eddy viscosity is then recovered via

ν t = ν f v1 , (8)
where the damping term f v1 is defined as

f v1 = χ 3 χ 3 + c 3 v1 with χ = ν ν . (9)
In the destruction term of Equation ( 7), the function f ω allows for a fast decay outside the boundary layer region, and is defined as

f ω (g) = g 1 + c 6 ω3 g 6 + c 6 ω3 (1/6) , g = r + c ω2 (r 6 -r) , r = ν SK 2 d 2 ω . ( 10
)
The distance-like parameter d ω defines a trigger which switches from this turbulence modeling in the boundary layer to a usual sub-grid scale viscosity in the rest of the simulation, and is defined as

d ω = min(d, C DES ∆x) ( 11
)
where d is the distance to the wall and ∆x is the local mesh size. Finally, the constants of the model are

c b1 = 0.1355 , c b2 = 0.622 , σ = 2 3 , K = 0.41 , c ω1 = c b1 K 2 + 1 + c b2 σ , c ω2 = 0.3 , c ω3 , c v1 = 7.1 . (12)

III. Application to a symmetric aerofoil case

To assess the capability of the turbulence modeling at solid boundaries to properly reproduce boundary layer dynamics at high Reynolds numbers, and therefore self-noise generation, a first benchmark case is studied. It consists in a symmetric NACA0012 aerofoil, at a temperature of 300K, a Mach number of 0.15 and a Reynolds number based on the chord c of the aerofoil of 6 million. The Reynolds number is defined as Re = u ∞ L/ν with u ∞ the velocity at infinity, L a characteristic length of the problem, and ν the kinematic viscosity.

This validation case has been extensively studied both numerically and experimentally [START_REF] Dennis | OVERFLOW Turbulence Modeling Resource Validation Results[END_REF], at different angle of attack (AoA), and is often used to assess the behaviour of a turbulence model [START_REF] Dennis | OVERFLOW Turbulence Modeling Resource Validation Results[END_REF].

A. Set-up

Based on the Mach number, the Reynolds number and the temperature of the benchmark, if we assume dry air, we can deduce the velocity at infinity u ∞ = 52.09m/s, the kinematic viscosity ν = 8.7.10 -6 m 2 /s and the pressure

P ∞ = 101364 Pa.
The mesh is composed of 6 levels of refinements, with a ratio of 2 between the mesh size of two adjacent refinement levels. Moreover, in this numerical approach, as the local CFL number is 1 in each grid cell, there is a direct link between the time step and the mesh size of each element. Therefore a multi time step multi mesh size method is used and only the finer elements are computed at each time step. The smallest mesh elements are located at the solid boundary, where the mesh size is 3.5.10 -4 m, which gives a local time step of 5.8.10 -7 s. A sketch of the refinement levels and their associated mesh size/time step used in this work is plotted in Figure 2.

The span of the aerfoil is set to 0.012 m, which correspond to a little over 1% of the chord. To limit blockage effects on the pressure due a restricted domain size, the computational domain extends to 20 chords in all direction except in the downstream direction where the domain extends to 40 chords, as plotted in Figure 2. The final mesh contains approximately 25 million cells. While this is a large number, only a fraction of these elements is computed at each time step (approximately 3 million).

B. Results

The computation is ran over 5 × 10 5 iterations, which corresponds to about 0.29 s of simulated time. This equates to the convection of a particle over 15 chords length. Each computation takes about 10h on 128 cores of Intel Nehalem X5560.

At no angle of attack, the evolution of the turbulent eddy viscosity ν t obtained from the Spalart-Allmaras turbulence model within the simulation is plotted in Figure 3a. Its thickness increases in the chord-wise direction of the aerofoil, which is coherent with the development of a boundary layer where more turbulence implies additional eddy viscosity. Moreover, it displays an expected symmetry between the upper and lower surface. In the Spalart-Allmaras model, the damping function f v1 in Equation 9 is defined such that ν (= ν χ) is linear in the viscous sub-layer, the buffer layer and the log-layer, therefore for y + (= u τ y/ν) up to a few hundreds. This behaviour is verified at multiple locations around the aerofoil, and is plotted for x/c = 0.3 on the case with no angle of attack in Figure 3b.

The pressure around the aerofoil, plotted in Figure 4a as the difference between the mean pressure p 0 and the pressure at infinity p ∞ , also displays an expected symmetry between the upper and lower surface, which is confirms in Figure 4b by the pressure coefficient C p defined as

C p = p -p ∞ 1 2 ρ ∞ u 2 ∞ . ( 13
)
The maximum pressure, located at the leading edge, is coherent with the physics, with a local pressure coefficient of 1 in the stagnation region. Similarly, as expected, the lowest pressure is observed on the upper and lower surfaces where the velocity experiences the maximum acceleration. Quantitatively and qualitatively, these results in Figure 4b are in very good agreement with the numerical prediction provided within the NASA benchmark [START_REF] Dennis | OVERFLOW Turbulence Modeling Resource Validation Results[END_REF].

IV. Application to a rod-aerofoil case

Interaction noise is an important aspect of turbomachine noise [START_REF] Khorrami | Towards full aircraft airframe noise prediction: lattice Boltzmann simulations[END_REF]. In this section, to evaluate this capability, a rod-aerofoil configuration is investigated. It is considered a representative benchmark [START_REF] Jiang | Numerical investigation on body-wake flow interaction over rod-airfoil configuration[END_REF] for numerical modeling of turbulent fluid motions interacting with airframe elements and is recommended by the Advisory Group for Aerospace Research and Development (AGARD). The rod creates a turbulent wake which is convected downstream and impinges the leading edge of the aerofoil. In this case, the boundary layer separation on the rod occurs upstream of the laminar to turbulent transition point which results in fully 3D turbulent structures in the wake. It forms a von Karman vortex street which creates an asymmetrical pressure distribution around the rod, producing periodic lift variations, therefore generating acoustic perturbations. Similarly, the unsteady lift on the aerofoil created by the incoming turbulence generates leading edge noise which also radiates to the far-field.

A. Set-up

A symmetric NACA0012 aerofoil with a chord c of 0.1 m is placed one chord-length downstream of a rod of diameter D = 0.1c, which corresponds to 0.01 m. The situation is depicted in Figure 5. The Reynolds number based on the diameter D of the rod is set to 4.8.10 4 , with an incoming velocity u ∞ = 72 m, a temperature T = 293 K and a mass density of the fluid ρ = 1.204 kg/m 3 . The speed of sound is a ∞ = 343 m/s giving a Mach number of approximately

M ≈ 0.2.
The mesh is composed of 6 levels of refinements, with a ratio of 2 between the mesh size of two adjacent refinement levels. The smallest mesh elements are located at the solid boundary, where the mesh size is 1.10 -4 m, which gives a local time step of 1.68.10 -7 s. This corresponds to a wall-distance in wall units, denoted x + , of about 20. A sketch of the refinement levels and their associated mesh size/time step used in this work is plotted in Figure 6.

B. Results

The computation is ran over 1.2 × 10 6 iterations, which corresponds to about 0.2 s of simulated time. This equates to the convection of a particle over 140 chords length. Each computation takes about 30h on 280 cores of Intel Xeon Broadwell E5-2680v4.

In Figure 7, a snapshot of the iso-surface of the Q-criterion colored by the streamwise velocity is plotted. It shows the boundary layer separation on the rod, creating a turbulent wake which visibly oscillates at a specific frequency, forming von Karman vortex streets. Then, the turbulent structures of the rod wake impinge the leading edge of the aerofoil, distorting the shape and sizes of these structures. At a point P1, located downstream of the rod, inside the rod wake (see Figure 5), the unsteady macroscopic values are recorded to perform a spectral analysis on the wake characteristics. The Power Spectral Density (PSD) on the streamwise velocity plotted in Figure 8a depicts a dominant oscillation at a Strouhal number * of S t = 0.192, which is in a very good agreement with experiments [START_REF] Jacob | A rod-airfoil experiment as a benchmark for broadband noise modeling[END_REF] and other numerical approaches [START_REF] Touil | Large-eddy simulation of rod-airfoil flow configuration in the subcritical turbulent regime by using the lattice Boltzmann method[END_REF]. Moreover, the high frequency part of the spectrum displays an expected Kolmogorov behaviour for the kinetic energy with a decay slope of -5/3.

The unsteady lift on the aerofoil surface generated by the wake interaction, produces an unsteady acoustic pressure which radiates to the far-field, as visible on the contour plot of the unsteady pressure in Figure 8b. From the Strouhal number of the vortex streets in the wake, the frequency of the acoustic perturbations are approximately 1380Hz, which results in a wavelength in the transverse direction of about 2.5c, which is in agreement with the unsteady pressure. The wavelength of these oscillations is large with respect to the chord of the aerofoil, therefore it radiates as a dipole. However, an interference pattern is visible in Figure 8b and originates from the unsteady lift on the rod, which also generates acoustical perturbations at the same frequency. These two acoustical sources being separated by a distance different than the wavelength, the two sources radiate with a shifted phase, creating interference patterns. Furthermore, in the far-field, the acoustic perturbations seem to travel unaffected by the change of mesh refinement and do not show any blocking effect due to a restricted domain size.

Finally, the mean pressure coefficients on the rod and on the aerofoil are plotted in Figure 9a and Figure 9b. As expected, they are perfectly symmetric. There is a good agreement with respect to the finite-volume Navier-Stokes LES (AVBP). Note that in this reference case, the LES is wall-resolved, therefore it does not rely on a turbulence model at the surface.

V. Application to a mobile turbofan engine case

With the objective of simulating complex turbofan applications, a more representative test-case which combines self-noise and interaction noise on mobile surfaces is investigated using a turbofan in rotation interacting with a grid-generated incoming turbulence.

A collaboration with Thomas Carolus of the University of Siegen is being set-up to perform experiments on this configuration. It will likely be based on a geometry already experimentally studied [START_REF] Carolus | A low pressure axial fan for benchmarking prediction methods for aerodynamic performance and sound[END_REF] but for self-noise, and will therefore be extended to consider interaction noise. A similar case has been studied experimentally [START_REF] Reese | Large eddy simulation of acoustical sources in a low pressure axial-flow fan encountering highly turbulent inflow[END_REF] and their set up is described in Figure 10.

This case is currently being set up in our solver, therefore there are no result yet. * Strouhal number S t =

3 f 3 c

 33 (x, t, c)dc and ρu(x, t) = ∫ f (x, t, c)dc.

Fig. 1

 1 Fig. 1 Representation of the particle velocities c α at each node of the D3Q19 lattice. The black node ( C ) points to the center of the lattice, the brown nodes ( F ) point to the faces of the neighboor cubes and the blue ( E ) and red nodes ( E ) point to the edges of the neighboor cubes. For clarity, only the velocities in the horizontal plane ( ) are displayed.

Fig. 2

 2 Fig. 2 Sketch of the mesh strategy of the symmetric aerofoil case. Each embedded colored-box corresponds to a different refinement levels, with the corresponding spatial and time step listed in the table below the sketch.

Fig. 3 (

 3 Fig. 3 (a) Contour of the turbulent eddy viscosity-like ν around the aerofoil obtained from the Spalart-Allmaras model (b) Evolution of ν at x/c = 0.3 as a function of y + , along with the linear relation ν = K yu τ ν imposed in the Spalart-Allmaras model from the viscous sub-layer to the log-layer via the function f v1 .

Fig. 4 ( 5 Fig. 5

 455 Fig. 4 (a) Contour of the difference between mean pressure p 0 and the pressure at infinity p inf (b) Pressure coefficient around the aerofoil.

Fig. 6

 6 Fig. 6 Sketch of the mesh strategy of the rod-aerofoil case. Each embedded colored-box corresponds to a different refinement levels, with the corresponding spatial and time step listed in the table below the sketch.

Fig. 8 (

 8 Fig. 8 (a) Strouhal number as a function of the power spectral density of the streamwise velocity at point P1 in the wake of the rod (b) Contour plot of the unsteady pressure perturbations in the simulation domain.

Fig. 9

 9 Fig. 9 Mean pressure coefficient (a) on the rod and (b) on the aerofoil.

VI. Content of the full paper

The full paper will include • a more detailed description of the numerical implementation of the DDES approach used in this work • regarding the NACA0012 aerofoil case, multiple angles of attack, including a separation case, will be investigated • a more comprehensive analysis of the rod-aerofoil case, with comparison with experiments and other numerical approaches of the separation bubble behind the rod, the RMS pressure coefficients of the rod/aerofoil, the velocity profiles at different locations, ... Also, it will include a comparison of the direct far-field against porous and solid Ffowcs Williams and Hawkings analogy • regarding the turbofan engine case, it will include the set-up, a detailed aerodynamic and aeroacoustic analysis, with comparisons with experiments.


Fran De Aquino
email: deaquino@elointernet.com

The Observed Magnitude of Cosmological Constant directly derived from the Quantization Theory of Gravity

Keywords: Cosmological constant, Cosmological constant problem, Quantization of Gravity, Quantization of Mass, Cosmology

Starting from the theoretical background developed for obtaining the Quantization of Gravity, in the paper Mathematical Foundations of the Relativistic Theory of Quantum Gravity, I show here that, it is possible to derive a theoretical value of the Cosmological Constant, extremely close to the observed value.

INTRODUCTION

The foremost challenge of Quantum Gravity Theory is to find the theoretical background, which would allow determining the observed positive value of Cosmological Constant Λ * , ( ) ( )

1 10 38 . 1 123 2 - × = Λ
observed planck l related to the accelerating expansion of our Universe [START_REF] Greene | The Hidden Reality: Parallel Universes and the Deep Laws of the Cosmos[END_REF][START_REF] Carroll | Spacetime and Geometry: An Introduction to General Relativity[END_REF]. The finding that the Universe is in accelerating expansion was first announced in 1998, and it results from astronomical observations of type IA supernovae carried out by Supernovae Cosmology Project and High-z Supernova Search Teams. The unexpected large discrepancy between theoretical results and observed result for Λ still exists [START_REF] Padmanabhan | [END_REF]4,5,6].

In this article I show that, starting from the theoretical background developed for the Quantization of Gravity, in the paper Mathematical Foundations of the Relativistic Theory of Quantum Gravity [START_REF] De Aquino | Mathematical Foundations of the Relativistic Theory of Quantum Gravity[END_REF] it is possible to derive a theoretical value of the Cosmological Constant, extremely close to the observed value (Eq. 1). * The Cosmological Constant is defined as

Λ = 8π(G/c )ρ = κρ 2 vac
vac , where κ is Einstein's constant. The true dimension of Λ is a length -2 . Thus, is a dimensionless number.

2 planck l Λ

THEORY

The Quantization of Gravity [START_REF] De Aquino | Mathematical Foundations of the Relativistic Theory of Quantum Gravity[END_REF] shown in the Mathematical Foundations of the Relativistic Theory of Quantum Gravity (See page 6) predicts the existence of an elementary quantum of matter (smallest indivisible particle of matter) with inertial mass, , given by (min given by (See [START_REF] De Aquino | Mathematical Foundations of the Relativistic Theory of Quantum Gravity[END_REF] page 34) ( )

0 i m ) ( ) ( )
3 min 0 i U U m M n = and ( ) ( ) 4 min 0 U i U CUF V m n = ρ where
is the total mass of particles (elementary quantum of matter) in the Universe, and is the volume of the Universe.

U M U V 2
Since the space, in the Universe, is expanding, then the space traversed by a photon expands behind it during the voyage since the Big-Bang. Thus, considering this fact, it is possible to show that the current distance to the most distant object we can see is evaluated in 46 billion light-years ( ) [9]). Thus,

m
3 3 4 U U R V π = is the volume of the expanding Universe ( ), while m R U 26 10 4 . 4 × = 3 3 4 H H R V π =
is the volume of the not expanding Universe. Consequently, we can write that ( )

5 U U U V M = ρ and ( ) 6 8 3 2 G H V M H U crit π ρ = =
From Eq. ( 5) and ( 6), we get ( )

7 8 3 3 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = = U H U U U R R G H V M π ρ
Now, starting from Eq. ( 3) and (2), we can write that ( )

8 8 3 2 max min 0 cd M m M n U i U U h π = =
By substitution of given by Eq.( 4) into Eq. ( 8), we obtain

U n ( ) ( ) 9 8 3 2 max min 0 cd M m V U i U CUF h π ρ = whence ( ) ( ) ( ) 10
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Λ U U i planck l CUF V M c m d G c G c G 3 2 1 h 43 42 1 ρ π
Therefore, Eq. ( 12) can be rewritten in the following form ( ) ( )

13 8 3 4 4 min 0 max 2 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = Λ U U i planck V M c m d G l Substitution of U U V M
given by Eq. ( 7) into Eq. ( 13), yields

  the maximum "diameter" that the Universe can reach [ See page 34).It also predicts that the space in the Universe is totally filled with particle of this type, forming a Continuous Universal Quantum Fluid (CUF), containing particles, with density U n CUF ρ , respectively

=

  Substituting the values in Eq.[START_REF] Carroll | The Cosmological Constant, See p.7[END_REF], we obtain the density of Cosmological Quantum Vacuum, i.e.,

  that theoretical value of the Cosmological Constant here obtained is extremely close to the observed value.