Maxime Philippe
email: maxime.philippe@innosea.fr

Bruno Borgarino
email: bruno.borgarino@innosea.fr

Panagiotis Kotronis
email: panagiotis.kotronis@ec-nantes.fr

Guillaume Ducrozet
email: guillaume.ducrozet@ec-nantes.fr

An integrated approach for the

published or not. The documents may come

INTRODUCTION

As the offshore wind industry expands in Europe, a wide variety of foundation designs are considered. Gravity based foundations (GBF) have the advantage of being competitive for shallow waters, and have a limited environmental impact during installation.

The design and certification processes for offshore wind turbine foundations require a high number of simulations to be performed (it can easily exceed 3000). In the case of steel foundations made of slender members such as monopiles, jackets or tripods, simulations by the Finite Element Method (FEM) can quickly be carried on, using beam elements. The wave loadings are computed by the Morison method. On the other hand, the engineering of concrete GBF has several challenges:

-The substructure itself is not slender, thus the Morison method is not appropriate to calculate wave loads. Plus, such substructures are installed in limited water depths, at which the linear wave theory is not valid for large wavelengths. -The use of solid, 3D finite elements to simulate reinforced concrete sections leads to high computational times. -The use of solid, 3D finite elements to simulate soil structure interaction can again lead to high computational times.

In the field of earthquake engineering, fast and accurate numerical tools have been recently developed, able to simulate the nonlinear behaviour of complex structures. Multifiber beam elements [START_REF] Kotronis | Simplified modelling strategies to simulate the dynamic behaviour of r/c walls[END_REF], [START_REF] Mazars | Using multifiber beams to account for shear and torsion. applications to concrete structural elements[END_REF], following the Euler-Bernoulli or the Timoshenko kinematic hypothesis, allow a simpler representation of the structure with a limited number of elements. Each section is made of different materials and each material follows a specific nonlinear constitutive law. This type of modelling strategy has been used to model the dynamic behaviour of reinforced concrete buildings submitted to earthquake loadings (see for example the 5-floor building simulated in [START_REF] Kotronis | Simplified modelling strategies to simulate the dynamic behaviour of r/c walls[END_REF], with only 24 beams and 37 fibers in each section). More recently, a 3-pier concrete viaduct has been modelled using 21 Timoshenko multifiber beam elements. Each section is made of 40 concrete fibers and 8 steel fibers [START_REF] Grange | The effects of soilstructure interaction on a reinforced concrete viaduct[END_REF]. For the dynamic simulations, concentrated (lumped) masses are introduced at the nodes.

When it comes to concrete structures under cyclic or earthquake loadings, suitable constitutive laws have to be considered. In the finite element calculations presented in [START_REF] Kotronis | Simplified modelling strategies to simulate the dynamic behaviour of r/c walls[END_REF] or in [START_REF] Grange | The effects of soilstructure interaction on a reinforced concrete viaduct[END_REF], the concrete and steel constitutive laws are chosen following [START_REF] Borderie | Phénomnes unilatéraux dans un matériau endommageable : modélisation et application à l'analyse des structures en béton[END_REF] and [START_REF] Menegotto | Method of analysis of cyclically loaded reinforced concrete plane frames including changes in geometry and non-elastic behaviour of elements under combined normal force and bending[END_REF] respectively. Concrete substructures are stiff, so most of the loads are transferred to the soil. Consequently soil-structure interactions must play a role in the foundations dynamics, and have to be taken into account in the analysis. The classical engineering approach to model these interactions is to use linear springs but their stiffness is difficult to calibrate a priori. On the other hand, a 3D finite element modelling to consider soil-structure interaction is time consuming and cannot be used for a large set of dynamic load cases. [START_REF] Grange | A macroelement for a circular foundation to simulate 3d soilstructure interaction[END_REF] developed a 2-nodes macro-element capable of representing the nonlinear soil structure interaction for a circular, superficial foundation. This macro-element is able to represent the following phenomena: elasticity, plasticity of the soil and uplift of the foundation. Some of the macro-element parameters can be calibrated using the soil properties and the foundation geometry. The other parameters do not vary much according to literature. It can be coupled with a multifiber finite element modelling for the upper structure, as shown in [START_REF] Grange | The effects of soilstructure interaction on a reinforced concrete viaduct[END_REF].

This paper aims at evaluating the applicability of these tools, to simulate the dynamic behaviour of a concrete GBF, with external forcing from the wind turbine, and wave loads computed with a nonlinear irregular wave model. First, a comparison is made for the case of a cantilevered substructure, with linear and nonlinear constitutive laws for the reinforced concrete. The objective is to demonstrate the ability of the model to capture nonlinear structural phenomena. Then, the nonlinear foundation macro-element is added to represent soil-structure interactions and to evaluate their effect on the GBF dynamics.

MODEL PROPERTIES Gravity based foundation properties

A general GBF design is needed to evaluate the applicability of the simulation tools. A simplified GBF design has been proposed, inspired by GBF substructures at the Thornton bank offshore wind farm [START_REF] Piere | Gravity base foundations for the thornton bank offshore wind farm[END_REF]. This substructure is by no mean a realistic, applicable design. The simplified GBF model is 45m-high and is situated in a water depth of 30 m. The foundation diameter is 20 m. The substructure is made of reinforced concrete and is ballasted (fully filled) with sand. Figure 1 shows the geometry of the GBF, and its main properties are summarized in Table 1

Wind turbine properties

The wind turbine used in this study is the reference turbine NREL 5MW [START_REF] Jonkman | Definition of a 5-MW Reference Wind Turbine for Offshore System Development[END_REF]. The hub-height is 95m and the tower height is 76.6m. This tower is the same as in the floating wind-turbine version of the NREL 5MW [START_REF] Jonkman | Definition of a floating system for phase IV of OC3[END_REF]. The connection between the substructure and the tower base occurs 15m above MSL (Mean Sea Level), at the top of the substructure.

Environmental conditions

The offshore wind turbine is subjected to turbulent wind and irregular waves loading. Wind and waves are aligned and propagate from the negative x-direction. Significant wave height of the considered sea state is 3m, and its peak period is 8s. The mean wind speed is 11.4 m.s -1 . The 600s time series of the irregular, 3D wind speed is based on a standard IEC Kaimal [10] turbulence model. Figure 2 shows the time series of the main component of the wind speed, in the x-direction. Wave elevation at the GBF location is represented on Figure 3.

METHODOLOGY Overview

The principle of the methodology presented in this paper is to integrate all relevant phenomena (aerodynamics, hydrodynamics, structural engineering, geotechnics) within the same simulation. Such a methodology leads to a more accurate combination of loads, thus reducing the sources of uncertainties. -FAST design code from NREL has been used to compute the wind turbine dynamics. -Hydrodynamic tools from LHEEA Lab. Ecole Centrale Nantes (ECN) have been used to compute hydrodynamic loads on the GBF. In-house nonlinear HOS incident wave model is used to compute wave kinematics. In-house diffraction/radiation code Aquaplus is used to compute linear diffraction loads. -Structural dynamics of the sub-structure is modelled with a nonlinear multifiber beam model from GeM Lab. Ecole Centrale Nantes. Soil-structure interaction is also computed with an in-house model from GeM Lab.

Numerical codes are described in further details in the following sections. The innovation of this study is the combination of these codes for offshore wind turbine applications.

Wind turbine loads and wave loads are computed in a first step, then applied to the substructure. Wind loads are applied at the tower-substructure connection. Wave loads are distributed over the nodes of finite element model of the substructure. In this paper no couplings between the wind and wave loads and the substructure are considered. The motions of the substructure are supposed small enough to consider it fixed. As a consequence:

-The wind turbine tower behaves as a cantilevered beam and the wind loads can be computed regardless the substructure; -No radiation hydrodynamic forces are applied to the substructure. The only wave loads are the Froude-Krylov excitation force and the diffraction force.

Several studies revealed that for flexible steel substructures such as monopiles, jackets or tripods the coupling between environmental loads and structural dynamics must be considered [START_REF] Jonkman | Offshore code comparison collaboration within iea wind task 23: Phase iv results regarding floating wind turbine modeling[END_REF][START_REF] Seidel | Integrated analysis of wind and wave loading for complex support structures of offshore wind turbines[END_REF][START_REF] Kaufer | Integrated analysis of the dynamics of offshore wind turbines with arbitrary support structures[END_REF]. To our knowledge, no results are available in the literature to determine if such coupling is mandatory for GBF. If necessary, coupling between the GBF and the turbine can be implemented in a future work, by using one of the following methods used for steel substructures:

-Sequential coupling [START_REF] Seidel | Integrated analysis of wind and wave loading for complex support structures of offshore wind turbines[END_REF], based on a linear condensation of the substructure finite element model. This method is easy to implement, but nonlinear structural and soil dynamics cannot be considered in the analysis. -Full coupling [START_REF] Kaufer | Integrated analysis of the dynamics of offshore wind turbines with arbitrary support structures[END_REF], based on the integration of the wind turbine representation in the finite element model as a macroelement. This method is more complex to manage, but is suitable for nonlinear calculations leading to a more accurate modelling of the system dynamics.

Hydrodynamic modelling

GBF are massive structures. Wave diffraction forces are significant and must be taken into account in the hydrodynamic calculations. Linear diffraction analysis supposes the linear Airy wave theory, which is only valid for small wave steepness and relatively deep water.

Henderson et al. [START_REF] Henderson | Hydrodynamic loading on offshore wind turbines[END_REF][START_REF] Henderson | Hydrodynamic loading on offshore wind turbine support structures[END_REF] studied the applicability of the hydrodynamic theory for the calculation of hydrodynamic loads on wind turbine substructures, in particular GBF. They concluded that it is necessary to evaluate the added mass coefficient of the support structure, and then to calculate the hydrodynamic loads based on the Froude-Krylov (pressure integration) method, using the added mass coefficient previously calculated in the diffraction analysis. The pressure integration method allows the use of a nonlinear wave theory.

Following these recommendations, a pressure integration method is used and combined with a nonlinear irregular wave theory to compute the Froude-Krylov loads. Diffraction loads are calculated following the linear potential theory.

The nonlinear wave model adopted is the High Order Spectral (HOS) wavemodel, developed at LHEEA Lab. The HOS method was initially proposed by West et al. [START_REF] West | A new numerical method for surface hydrodynamics[END_REF] and Dommermuth and Yue [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. The development of this method at LHEEA Lab. lead to a numerical wave tank (including wave generation, absorption, wall reflexions) [START_REF] Ducrozet | A modified high-order spectral method for wavemaker modeling in a numerical wave tank[END_REF] and an oceanic wave model [START_REF] Ducrozet | 3-d hos simulations of extreme waves in open seas[END_REF], which can be coupled to sea-keeping models. The methodology for coupling the linear potential sea-keeping model and the HOS model [START_REF] Gilloteaux | A non-linear potential model to predict large-amplitude-motions: Application to the searev wave energy converter[END_REF] has been used in the present study.

Froude-Krylov loads are computed by integrating the dynamic pressure, calculated by the HOS model, on the instantaneous wetted surface. Diffraction loads are calculated using a combination of the linear potential method and the wave elevation calculated with the HOS model (see eg. Equation 3). Hydrostatic loads are computed as a mean buoyancy force applied on the nodes. One limit of the present methodology is that it does not take into account breaking wave loads.

Hydrodynamic loads are computed by the Boundary Element Method. Velocity potential, hydrodynamic pressure and the resulting forces are expressed at the center of each panel of a surface mesh of the substructure (see figure 6). These distributed forces have to be at the nodes of the finite element model. For the Froude-Krylov force and the diffraction force, the following method is applied:

-Compute the hydrodynamic force at each panel; -For each node i of the finite element mesh, select the N i panels of the BEM (Boundary Element Method) mesh whose centers z k satisfy:

(z i-1 -z i )/2 < z k ≤ (z i -z i+1 )/2 (1)
-Sum the contributions of the N i panels at each node i:

F i = N i ∑ k=1 F k (2)
The wave loads are computed in the time domain, for each node i and degree of freedom q, using the following equation:

F q ex i (t) = N i ∑ k=1 p dyn n q k S k Froude-Krylov + +∞ -∞ K q dif i (t -τ)η(t)dτ diffraction ( 3
)
p dyn is the dynamic pressure due to the incident wave field, and η is the surface elevation at the substructure center. p dyn and η are calculated by the nonlinear HOS wave model. n i is the i-th component of the generalized normal vector of the panel k, S k is the surface of this panel. K q dif i is the diffraction impulse response. This response is obtained based on RAOs (Response Amplitude Operator) from the linear, frequency-domain hydrodynamic software Aquaplus [START_REF] Delhommeau | Seakeeping codes aquadyn and aquaplus[END_REF], using equation 4:

K q dif i (t) = 1 2π +∞ -∞ f q dif i (ω)e iωt dω (4)
In Equation 4, f q dif i is the RAO of the diffraction force at node i, computed from Aquaplus results on each panel using Equation 2. Figure 4 shows these values computed at nodes 2 to 5. It can be seen that forces are highest at node 5, because it is loaded on a large vertical span (from the free surface to the midpoint between nodes 4 and 5, see Figure 1). Nodes 3 and 4 come second because the substructure is larger at their level, thus diffraction loads are more important. Node 4 is the less loaded, probably because it is situated at mid-depth, where the substructure is only moderately large.

Figure 5 shows the final time domain results for the Froude-Krylov and the diffraction forces in the surge direction.

Structural modelling

A finite element model using multifiber beams and concentrated masses is used to reproduce the nonlinear structural response of the GBF. The structure is discretised using 6 Bernoulli multifiber beams represented on Figure 1. 120 concrete fibers and 216 steel fibers are used in each section. Steel fibers represent the reinforcement bars at their actual position (regularly positioned in 2 concentric circles). Hydrodynamic loads and wind loads are computed independently and applied at the nodes. The weight of concrete, steel and ballast are also applied on the nodes. Calculations are made with FEDEASLab, a finite element matlab MATLAB toolbox [START_REF] Fc Filippou | FEDEASLab getting started guideand simulations examples[END_REF].

The concrete constitutive model is based on damage me- chanics and is able to reproduce the nonlinear behaviour of the material for cyclic or dynamic loadings [START_REF] Borderie | Phénomnes unilatéraux dans un matériau endommageable : modélisation et application à l'analyse des structures en béton[END_REF]. The plasticity model adopted for steel is a modified version of the classical Menegotto-Pinto model [START_REF] Menegotto | Method of analysis of cyclically loaded reinforced concrete plane frames including changes in geometry and non-elastic behaviour of elements under combined normal force and bending[END_REF]. Material data for concrete and steel are listed in table 2. E is the Young modulus, ν is the Poisson ratio. Y i are the associated variables to the damage variable (i = 1 for tension or 2 for compression) and Y 0i are the initial elastic thresholds. A i , B i and β i are material constants. σ f is the stress for crack closure. f y and ε sh are the steel stress and strain at yielding and f su and ε su at rupture.

Soil-structure interaction modelling

A macro-element model is introduced at the base of the finite element model to take into account soil structure interaction. The macro-element is able to simulate the 3D behaviour of a rigid shallow foundation, submitted to cyclic loadings. The general formulation of the macro-element is detailed in [START_REF] Grange | A macroelement for a circular foundation to simulate 3d soilstructure interaction[END_REF] for the plasticity model, [START_REF] Grange | A macro-element for a shallow foundation to simulate soilstructure interaction considering uplift[END_REF] for the uplift mechanism, [START_REF] Grange | A macroelement to simulate 3d soil-structure interaction considering plasticity and uplift[END_REF] for the coupling of the two mechanics and [START_REF] Grange | A macroelement to simulate dynamic soil-structure interaction[END_REF] for the dynamic behaviour.

In the following, a class B soil according to the Eurocode 8 classification is considered at the base of the GBF. The adopted macro-element parameters are detailed in [START_REF] Grange | The effects of soilstructure interaction on a reinforced concrete viaduct[END_REF], see table 3. As mentioned before, the aim of the article is to show the ability of the macro-element to model the soil structure interaction for a GBF following an integrated approach. For a more precise study, the macro-element parameters have to be identified for a marine soil.

Wind turbine modelling

In order to accurately assess the effect of wind loads on the structure, the aeroelastic software FAST [START_REF] Jonkman | FAST User Guide[END_REF] is used. FAST has the capability to fully model the dynamics of a wind turbine, including elasticity, aerodynamics and control in the time domain. As a consequence, the effects of cyclic loadings due to turbulent wind and of aerodynamic damping can be modelled, which removes possible sources of uncertainties and over-conservatism. Loads computed at the bottom of the turbine tower using FAST are then applied in the finite element model of the GBF substructure. Figure 7 shows the wind turbines loads, computed with FAST, and applied at the top of the substructure in the xdirection. The NREL 5 MW wind turbine is controlled by a pitch controller. A mean wind speed of 11.4 m/s has been chosen. This case is critical as the rotor thrust is maximal for this wind speed [START_REF] Jonkman | Definition of a 5-MW Reference Wind Turbine for Offshore System Development[END_REF].

RESULTS

Effect of nonlinear constitutive laws for steel and concrete

In this section the results of two different finite element calculations are compared, considering linear and nonlinear constitutive laws for steel and concrete of the reinforced concrete structure. Soil structure interaction is not taken into account as the GBF is supposed fixed at its base. The aim is not to evaluate the design of the foundation, which is not realistic, but to evaluate the applicability of the proposed integrated methodology to reproduce nonlinear phenomena such as decrease in the structural stiffness due to damage and permanent strains.

Top horizontal displacements: Figure 8 shows the GBF top horizontal displacements considering a linear 8a and a nonlinear structural modelling 8b. One can observe that the mean deflection is found more important for the nonlinear structure. In other words, for this level of loadings and without considering soil structure interaction structural stiffness has decreased.

Concrete axial stresses:

This decrease of the structural stiffness implies not also an increase of the GBF top displacements but also a decrease of the applied concrete stresses in the damaged sections. Figure 9 shows the axial concrete stresses in a section at the level of node 5 ; figure 8a for linear constitutive laws, and figure 8b for nonlinear constitutive laws. We can observe a decrease in axial concrete stresses in the non linear simulation. The section is damaged, cracking appears. Load is transmitted to the steel reinforcement bars that should be able to deal with more important axial stresses. 4 that stresses in reinforcement bars actually increase when calculated using the non linear structural model. Nevertheless steel remains in the elastic domain. Figure 10 shows the axial stresses in a steel bar in a damaged section at the level of node 5, calculated with the non linear structural model.

Moment at the base:

Even if nonlinear phenomena appear on the top of the substructure, it does not seem to have significant influence on the way forces are transmitted on the foundation. This is probably due to the fact that nonlinearities are limited (reinforcement remains elastic).

The power spectral density of the stochastic time series of the moment is shown in Figure 11. Such representation allows a first identification of the phenomena which influence the response. The main peak at 0.12 Hz corresponds to wave loads (the peak period being 8 s). Power at lower frequencies are due to wind loads. For this loading case, waves slightly dominate the moment response of the substructure. At 11.4 m/s, the rotational speed of the turbine is 12 rpm. The 1P turbine excitation frequency is thus 0.2 Hz, and the 3P frequency is 0.6 Hz (1p equals the rotational frequency, and 3p equals 3 times the rotational frequency). The 3P frequency is clearly observable on the graph. It is also interesting to note a peak around 0.4 Hz, probably corresponding to the first natural frequency of the turbine tower. This shows that turbine loads are properly transmitted from the turbine to the substructure base.

Damage: Damage due to traction is only observed at the levels of nodes 5 and 6 (the rest of the structure remains in the linear regime), see Figure 12. No damage due to compression is found. One half of the section is damaged as the loading is important only in one direction (not a "true" symmetric load with significant maximum values at both directions).

Effects of soil-structure interaction

To assess the ability of the integrated numerical strategy to reproduce the effects of soil structure interaction, we compared hereafter the numerical results of the GFB considering or not the macro-element. For both calculations, the nonlinear constitutive laws for concrete and steel are used.

Table 5 shows some statistics of the two simulations. It is found that the horizontal top displacements of the GBF top are significantly increased, in terms of mean and maximum displacements. The results are however similar in terms of moment at the foundation base and stresses in the steel bars. As in [START_REF] Grange | The effects of soilstructure interaction on a reinforced concrete viaduct[END_REF], the GBF on a B type soil seems to behave more as a fixed structure with the development of plastic hinges in the substructure. Results could be significant different if a less rigid type soil (C) was considered. Results presented here have to be confirmed with macro-element parameters calculated for a marine soil.

CONCLUSION AND FUTURE WORK

This paper shows the interest of an integrated modelling strategy combining multifiber finite element beams with nonlinear constitutive laws, nonlinear soil structure interaction (macroelement) and nonlinear wave modelling for the structural dynamic modelling of concrete gravity base foundations. The model has been able to efficiently simulate the behaviour of a GBF foundation under nonlinear irregular wave forcing and loads transmitted by the turbine. The model reproduces complex nonlinear phenomena such as decrease in the structural stiffness due to damage, permanent strains, and the amplification of the horizontal top displacement due to soil structure interaction.

Even if each numerical tool used in this study has already been validated independently, it will now be necessary to vali- date the combination of these tools for offshore wind turbine applications. Results have to be validated using more advanced and time consuming numerical models, including 3D FEM models for soil-structure interaction and CFD models for fluid-structure interactions. Comparison with experimental data would also be of great interest. Once validated, the model will be applied on a variety of design load cases including fatigue limit states and extreme limit states. This work will now continue in the following ways :

-A more detailed analysis of the structural dynamics of the GBF will be performed including comparison of multiple design load cases in order to identify the most severe ones for the concrete structure. The influence of wind loads and wave loads on the dynamics of the structure will be quantified and design guidelines will be provided. -A more refined mesh will be used for the nonlinear structural analysis in order to identify and to localize with more precision the existence and the position of the plastic hinges -Parametric studies will be performed to identify the influence of parameters such as wind speed, wave conditions and soil stiffness. Results will be compared to reference results on other type of substructure (such as the OC3 project [START_REF] Jonkman | Offshore code comparison collaboration within iea wind task 23: Phase iv results regarding floating wind turbine modeling[END_REF]), to assess the need of a coupling between the structure model and the wind turbine model. -The effect of soil-structure interaction will be assessed through comparison considering different types of soil.

FIGURE 1 :

 1 FIGURE 1: GBF Geometry. Measures are in mm. Blue dots represent the nodes of the multifiber finite element model; blue lines represent the finite elements. See Section "Structural modelling"

FIGURE 2 :

 2 FIGURE 2: Main component of the 3D turbulent wind field at hub-height

FIGURE 3 :

 3 FIGURE 3: Wave elevation at GBF location

FIGURE 4 :

 4 FIGURE 4: Diffraction force RAOs in surge, computed at FEM nodes (see Figure 1)

FIGURE 5 :

 5 FIGURE 5: Wave loads at node 2 in the surge direction. Red, solid line: Froude-Krylov excitation force. Blue, dotted line: diffraction force.

FIGURE 6 :

 6 FIGURE 6: Froude-Krylov loads on GBF

  θ θ = 1299456.1 MNm/rad q max = 2

FIGURE 7 :

 7 FIGURE 7: Wind turbine loads applied at foundation top in x direction

FIGURE 8 :

 8 FIGURE 8: Comparison of the GBF top horizontal displacements for the linear (fig. 8a) and the nonlinear (fig. 8b) structural analysis

FIGURE 9 :

 9 FIGURE 9: Comparison of concrete axial stresses in a section at the level of node 5 for the linear (fig. 8a) and the nonlinear (fig. 8b) structural analysis

FIGURE 10 :

 10 FIGURE 10: Stresses in a steel bar in a damage section at the level of node 5 (non linear structural analysis)

FIGURE 11 :

 11 FIGURE 11: Power spectral density on the moment at the substructure base (non linear structural analysis)

FIGURE 12 :

 12 FIGURE 12: Damage due to traction in a section at the level of node 5 (non linear structural analysis)

TABLE 1 :

 1 . Main properties of the GBF

	Base diameter	20 m
	Diameter at connection with turbine	6.5 m
	Foundation height	45 m
	Foundation draft	30 m
	Wall thickness	0.5 m
	Base height	1.9 m
	Number of steel bars in a section	216
	Steel bars diameter	4 cm
	Concrete mass	2653 t
	Steel mass	327 t
	Ballast mass	5445 t

TABLE 2 :

 2 Material data for concrete and steel

	Concrete parameters Steel parameters
	E	20 GPa	E	200 GPa
	nu	0.2	f y	460 MPa
	Y 01	380 Pa	f su	710 MPa
	Y 02	0.9 GPa	ε sh	0.0023
	A 1	6.00E-03	ε su	0.11
	A 2	5.00E-06
	B 1	1.0
	B 2	1.6
	β 1	1 MPa
	β 2	-40 MPa
	σ f	3.5 MPa

TABLE 3 :

 3 Macro-element parameters

TABLE 4 :

 4 Comparison of statistical results for the linear and the nonlinear structural analysis

TABLE 5 :

 5 Comparison of statistical results considering or not soil structure interaction (NL: NonLinear structural analysis and cantilevered foundation ; NL+SSI: Non Linear structural analysis and macro-element for Soil Structure Interaction)

ACKNOWLEDGMENT

The authors would like to acknowledge Kévin GUIOT from INNOSEA for his help in GBF geometry definition, and Rocio ISORNA from INNOSEA for her help in the structural analysis.


I Ribet-Mohamed

V Arounassalame

J Nghiem

M Caes

M Guénin

L Höglund

E Costard

P Christol

Residual fixed pattern noise and random telegraph signal noise of a MWIR T2SL focal plane array

Keywords: stability over time, FPN, RFPN, correctability, T2SL, MWIR, FPA

Stability over time has recently become a figure of merit of major importance to compare the performances of infrared focal plane arrays (FPA) of different technologies. Indeed, this parameter dictates how often the calibration of operational electro-optical systems has to be done, and thus reflects the availability of the system during an operational mission.

The stability over time is generally estimated through fixed pattern noise (FPN) and residual fixed pattern noise (RFPN) measurements after a two-point correction. However, each laboratory or industrial has its own protocols and criteria, such that published results cannot be easily compared. Recent studies also showed that random telegraph signal (RTS) noise, which leads to flickering pixels, can strongly affect the image quality, so the question arises as to wether these RTS pixels have an effect on RFPN.

In this paper, we describe our experimental protocol to evaluate the stability over time of an FPA and to count up / classify flickering pixels. We then present the results obtained on a T2SL MWIR Integrated Detector Dewar Cooler Assembly (IDDCA) provided by IRnova. Our measurements show that the stability over time of the T2SL MWIR IDDCA are excellent: first, in terms of FPN/RFPN; then, in terms of RTS noise with only a few blinking pixels. We also show that the RTS pixels having an effect on the RFPN are fully detected by the algorithm used to rule out defective pixels before calculating RFPN.

Introduction

Type-II Superlattice (T2SL) detectors are a recent commercial infrared technology designed to meet high-performance imaging applications [START_REF] Sood | Design and Development of Two-Dimensional Strained Layer Superlattice (SLS) Detector Arrays for IR Applications, Two-dimensional Materials for Photodetector Pramoda Kumar Nayak[END_REF][START_REF] Rogalski | InAs/GaSb type-II superlattice infrared detectors: three decades of development[END_REF][START_REF] Chevallier | Dark current reduction in microjunction-based compound electron barrier type-II InAs/InAs1-xSbx superlattice-based long-wavelength infrared photodetectors[END_REF][START_REF] Ting | Antimonide type-II superlattice barrier infrared detectors[END_REF][START_REF] Klipstein | Long Wave Infrared Type II Superlattice Focal Plane Array Detector[END_REF][START_REF] Höglund | Manufacturability of type-II InAs/GaSb superlattice detectors for infrared imaging[END_REF][START_REF] Durlin | InAs/InAsSb superlattice structure tailored for detection of the full midwave infrared spectral domain[END_REF][START_REF] Daumer | Photodetector development at Fraunhofer IAF: From LWIR to SWIR operating from cryogenic close to room temperature[END_REF][START_REF] Kazemi | Mid-wavelength infrared unipolar nBp superlattice photodetector[END_REF][START_REF] Tan | The 640 × 512 LWIR type-II superlattice detectors operating at 110 K[END_REF], such that electro-optical characterisation efforts are now focused on system-oriented merit functions, for example MTF (Modulation Transfer Function) or stability over time. MTF describes how well a detector can reproduce spatial frequencies [START_REF] Nghiem | MTF measurements of a type-II superlattice infrared focal plane array sealed in a cryocooler[END_REF], whereas stability over time dictates how often the calibration of operational electro-optical systems has to be done, and thus reflects the availability of the system during an operational mission.

The stability over time is generally estimated through fixed pattern noise (FPN) and residual fixed pattern noise (RFPN) measurements, the latter being evaluated after a linear two-point calibration (TPC). However, each laboratory or industrial has its own protocol and criteria to exclude some pixels considered as defective, such that published results cannot be easily compared. Moreover, one can wonder wether FPN and RFPN are affected by pixels with random telegraph signal (RTS) noise. The signal delivered by these pixels can oscillate between two (or more) levels at certain times, creating a particularly harmful blinking effect for image quality. Since the blinking is not permanent, detecting these pixels is a very difficult task.

In this paper, we describe our experimental protocol to evaluate the temporal stability of a FPA. We then present the results of long-term measurement campaigns realized on an InAs/GaSb T2SL focal plane array. We chose to carry out our measurements on a detector integrated in an operational packaging (IDDCA, which means Integrated Detector Dewar Cooler Assembly) in order to remain as close as possible to the real conditions of use. More precisely, our measurements were realized on a T2SL MWIR 320x256 pixels IDDCA provided by IRnova, on which we have already published very encouraging first stability over time measurements [START_REF] Nghiem | MTF and FPN measurements to evaluate midwave infrared T2SL focal plane arrays[END_REF].

The first paragraph of this paper is dedicated to residual fixed pattern noise (RFPN) measurements. After recalling the principle of this measurement, we present our experimental setup and the data processing. Results obtained during a long-term measurement campaign are then reported and commented. The second paragraph focuses on random telegraph signal (RTS) noise, and presents the experimental setup, the data processing and the experimental results.

Residual fixed pattern noise measurements

The stability over time is generally estimated through fixed pattern noise (FPN). The FPN corresponds to the spatial fluctuations of the signal delivered by the pixels of a focal plane array (FPA) when it receives a homogeneous incident power. The fixed pattern noise can originate in the pixel itself (for example if there are differences in cut-off wavelengths or variations in active area thickness through the focal plane array) or in the read-out integrated circuit (for example if the output amplifiers have slightly different characteristics (gain, offset)). Non-uniformity correction (NUC) is generally realized with a linear two-point calibration [START_REF] Perry | Linear theory of nonuniformity correction in infrared staring sensors[END_REF]. To implement it, two uniform backgrounds are presented to the FPA (generated for example with two extended blackbodies at different temperatures). Assuming that the signal delivered by each pixel varies linearly with the incident power, one can calculate, for each pixel, a gain correction coefficient G i,j and an offset correction coefficient O i,j , accounting for multiplicative and additive non-uniformities, respectively. They are defined as :

        2 , 1 , 2 , 1 , ,     j i j i j i j i j i S S S S G    (1)
    1 , , 1 , ,   j i j i j i j i S G S O    (2)
where S i,j ( 1 ) (resp. S i,j ( 2 )) is the signal delivered by the pixel (i,j) for an incident power equal to  1 (resp.  2 ) and after removing the ROIC offset ; <S i,j ( 1 )> (resp. <S i,j ( 2 )>) is the spatially averaged signal delivered by the FPA (after removing the ROIC offset) for an incident power equal to  1 (resp.  2 ).

After 2 points correction, the signal delivered by the pixel (i,j) becomes:

    j i j i j i j i G O S S , , , , '     (3)
When applying the two-point correction that has just been calculated, the corrected image is perfectly uniform at an incident power equal to  1 or  2 .

The two-point calibration has the advantage of being easy to implement and very effective at the time it is calculated. But the question is how long this calibration remains valid. To answer this question, the fixed spatial noise that persists after the two-point correction -called residual fixed spatial noise (RFPN) -is evaluated:

            i j j i j i S S N RFPN 2 , , 1    (4)
where N is the total number of pixels of the FPA.

To evaluate the performance of a system, temporal noise and residual fixed pattern noise are generally summed quadratically [START_REF] Mooney | Responsivity nonuniformity limited performance of infrared staring cameras[END_REF], the former being defined as the temporal fluctuations of the signal delivered by a pixel. Consequently, the influence of spatial noise decreases rapidly as soon as the RFPN is lower than the temporal noise; therefore, one can consider that the calibration is valid as long as the RFPN remains lower than the temporal noise, even if stricter criteria have been proposed for particular conditions [START_REF] Mooney | Characterizing IR FPA nonuniformity and IR camera spatial noise[END_REF] (with a human observer in the loop).

Experimental setup

The measuring bench relies on an extended area blackbody which is placed in front of the IDDCA. Scanning the blackbody temperature allows to scan the incident power on the focal plane array. For each blackbody temperature explored, a 256-image cube is recorded at an integration time of 4ms, and another one at null integration time to measure the read-out integrated circuit (ROIC) offset.

Data processing

The analysis of the measurements consists of three main stages:  the detection of defective pixels, and once these have been set aside,  the application of the two-point correction  the calculation of the residual fixed pattern noise In order to highlight the importance of the criteria used to define a defective pixel, we compared several algorithms. The first one, which will be further referred as the reference algorithm, considers the average image recorded for a 50% well fill and corrects it for the dome effect due to variation in the geometrical throughput between the center and the corner of the focal plane array. It then sets aside pixels with continuous level, noise (defined as the standard deviation of the signal of one pixel over 256 images), responsivity or noise equivalent temperature difference (NETD) too far from the average values. The exact criteria are summarized in Table 1, together with the number of pixels falling into each category. The total number of defective pixels detected with this algorithm is 41, which corresponds to an operability of 99.95%. As can be seen in Figure 2, the spatial distribution of the defective pixels is relatively homogeneous, with only 3 clusters with 2x1 pixels, and no bigger cluster.

Criteria

Number of pixels Continuous level

j i j i j i S S S , , , % 30    31 Noise j i j i j i , , , % 50       25 Responsivity j i j i j i R R R , , , % 30    34 NETD j i j i j i NETD NETD NETD , , , % 100    37 Total 41
Table 1 -Criteria used to classify one pixel as a defective one (reference algorithm). S i,j is the signal delivered by the pixel, < S i,j> is the spatially averaged signal,  i,j is the noise of the pixel, < i,j> is the spatially averaged noise, R i,j is the responsivity of the pixel, < R i,j> is the spatially averaged responsivity, NETD i,j is the NETD of the pixel and <NETD i,j > the spatially averaged NETD. The total number of defective pixels is lower than the sum of the number of pixels sorted as defective with only one criterion, since most of the defective pixels fall into more than one category.

The second algorithm, which we will call advanced algorithm, is based on the same principle, except that a two-point correction (calculated from images at 33% and 66% well fill) is applied after the dome effect is corrected and before the defective pixels are detected. The exact criteria are summarized in Table 2, together with the number of pixels falling into each category. The total number of defective pixels detected with this algorithm is 49, which corresponds to an operability of 99.94%. As can be seen on the defective pixels map presented in Figure 2, the spatial repartition of defective pixels is once again quite homogeneous.

Criteria

Number of pixels Continuous level Once the defective pixels have been set aside, we apply the gain and offset corrections (see equation ( 3)) and then calculate the residual fixed pattern noise thanks to equation (4). Results are presented in the following paragraph.

j i j i j i S S , , , 10     43 Noise j i j i j i , , , % 50       25 Responsivity j i j i j i R R R , , , % 30    34 NETD j i j i j i NETD NETD NETD , , , % 100    37 Total 49

Results and discussion

Figure 3 presents the ratio between the residual fixed pattern noise and the temporal noise (TN) as a function of the well fill factor. The left-hand side (resp. right-hand side) graph corresponds to the case where defective pixels are detected with the reference algorithm (resp. advanced algorithm). In both cases, the gain and offset corrections were calculated on the first day of the campaign. The RFPN/TN results obtained that day exhibit the expected W-shape, reaching zero for the two values of the power used to calculate the gain and offset coefficients. The RFPN/TN measurements obtained on the following days are all higher than those obtained on the first day, but if we trace their temporal evolution (see Figure 4), it appears that the RFPN is remarkably stable over 6 weeks. More precisely, the ratio between the RFPN and the temporal noise is typically equal to 1.1 for the reference algorithm and equal to 0.83 for the advanced algorithm. These results clearly show that the temporal stability of the T2SL MWIR FPA is excellent. They also stress the importance of the data processing phase and the criteria used to rule out defective pixels, since removing 8 more pixels decreases the RFPN/TN plateau of 25%. It appeared that these 8 pixels presented a signal close to the average for 50% well fill (and were therefore not detected with the reference algorithm), but that they were strongly non-linear and impossible to correct with a TPC calibration. That's why the second algorithm detected them.

Random telegraph signal noise study

Pixels are considered as exhibiting random telegraph signal (RTS) noise when their signal oscillates between at least two levels, while the received power and operating conditions do not change. RTS noise can be particularly detrimental to image quality since the resulting blinking effect cannot be suppressed by 2-point correction. It is thus considered as a major challenge that different technologies will have to face [START_REF] D'souza | Noise attributes of LWIR HDVIP HgCdTe detectors[END_REF][START_REF] Brunner | Improvement of RTS noise in HgCdTe MWIR detectors[END_REF][START_REF] Pogány | Random telegraph signal noise instabilities in lattice-mismatched InGaAs/InP photodiodes[END_REF], especially in addressing high operating (HOT) and size, weight and power (SWaP) requirements.

Experimental setup and data processing

Studying RTS noise requires long temporal acquisitions under fixed conditions (integration time, blackbody temperature…). We determined that recording cubes of 5000 images was a good compromise between RTS pixels detection efficiency and data processing complexity. The integration time was set to 4ms. In order to study the effect of RTS pixels on RFPN, we needed simultaneous RFPN measurements so we scanned the blackbody temperature between 10°C and 58°C by steps of 3°C, and each time recorded a 5000 images data cube. Such RFPN/RTS measurements were repeated 24 times, representing a cumulated acquisition time of more than 9 hours. Some of the measurements were made within the same cooling, others after allowing the detector to return to room temperature.

To detect RTS pixels, we used an algorithm which starts by isolating those with excess temporal noise. Among these candidates, it then fits the histogram with one or several Gaussian functions: if two or more Gaussian functions are required, the pixel is considered as affected by RTS noise (see example Figure 5).

Results and discussion

Our measurements show that the number of RTS pixels of the MWIR T2SL 320x256 pixels FPA is very low. It varies between 0 and 10 pixels detected by acquisition of 5000 images (with a mean value of 6 RTS pixel per cube (acquired in 80 seconds). Table 3 presents a classification of the 52 RTS pixels detected on 192 cubes (for each RFPN measurement, only cubes corresponding to blackbody temperatures between 10°C and 31°C were considered, to ensure that RTS pixels would not be hidden by an increasing photon noise). Only 19 pixels have a true 2-level behaviour (see Figure 5 a andb). 15 pixels exhibit spikes (see Figure 5 c). To finish with, 18 pixels exhibited low frequency temporal fluctuations that were considered as RTS noise (see Figure 5 d). We also plotted in Figure 6 the instantaneous and the cumulative number of RTS pixels detected as a function of the measurement number. The cumulative number of RTS pixels increases more rapidly at the beginning of the curve and less rapidly for measurement numbers > 125, which could be compatible with a stabilisation of the cumulative number of RTS pixels but the stabilisation regime is still not reached after 192 measurements (~4.5 hours). Figure 6 also shows that the increase of the total number of RTS pixels is not correlated with new cooling cycles (materialized by arrows). Finally, we wondered whether the RTS pixels could affect the RFPN. From the data cubes acquired at different blackbody temperatures, we extracted RFPN measurements using the protocol described above. We first ruled out the 53 pixels considered as defective by our "advanced" algorithm (including 10 RTS pixels) and computed the RFPN/TN ratio as a function of the well fill, see Figure 7). We then also ruled out the 6 other RTS pixels which were not detected as defective, and computed the RFPN/TN ratio. Both curves superimpose, which proves that the RTS pixels which affect the RFPN were already detected by our "advanced" algorithm : there is no need to improve it.

Conclusion

In this paper, we described our experimental protocol to evaluate the stability over time of an FPA and to count up / classify pixels with random telegraph signal noise. We presented the results obtained for a T2SL MWIR 320x256 pixels IDDCA provided by IRnova. The stability over time is excellent over seven weeks. We stressed out the importance of the criteria used to rule out defective pixels. The RFPN/TN plateau can thus be decreased from 1.1 (with reference defective pixels detection and basic two-point correction (TPC)) to 0.83 (with advanced defective pixel detection and basic TPC).

Our measurements also allowed us to count up and classify pixels exhibiting a random telegraph signal (RTS) noise. The total number of 2-level RTS pixels is 19, for a 4.55 hours cumulative acquisition. 15 out of these 19 2-level RTS pixels flicker less than 10% of the time. We also showed that among the RTS pixels, those affecting the RFPN were already detected as defective by our "advanced" algorithm.

Figure 1 -

 1 Figure 1 -T2SL MWIR IDDCA from IRnova.

Figure 2 -

 2 Figure 2 -Defective pixels map detected with the reference algorithm (black dots) and by the advanced algorithm (black dots + red circles). 41 (resp 49) pixels are set aside with the reference (resp advanced) algorithm, which represents an operability of 99.95 (resp 99.94%).

Figure 3 -Figure 4 -

 34 Figure 3 -Ratio between residual fixed pattern noise (RFPN) and temporal noise (TN) as a function of well fill, for measurements realized during 7 weeks. Left : 41 defective pixels (identified with the reference algorithm) were set aside ; Right : 49 defective pixels (identified with the improved algorithm) were set aside. The gain and offset coefficient used in the two-point corrections are those calculated on the first day of the campaign.

Figure 5 -

 5 Figure 5 -Left : signal delivered by a pixel affected by a 2-level RTS noise as a function of time ; Right : associated histogram, showing two distinct populations.

Table 3 -Table 4 -

 34 Classification of pixels detected as affected by RTS noise. Signal as a function of time (one frame is taken each 16.7ms) for pixels detected as RTS : a) and b) exhibit a clear 2 level behavior, whereas c) exhibits spikes and d) exhibits low frequency fluctuations. Let's now focus on the 19 pixels exhibiting a clear 2 level-behavior. 15 pixels out of them flicker less than 10% of the time (for example, pixel n°28430, shown in Figure 5 a, only flickers 2% of the time), 3 pixels flickers between 15% and 20% of the time and one pixel (n°5409, see Figure 5 b) flickers 90% of the time.

Figure 6 -

 6 Figure 6 -Instantaneous and cumulative number of RTS pixels detected as a function of the measurement number. Red arrows indicate new cooling cycles.

Figure 7 -

 7 Figure 7 -Ratio between residual fixed pattern noise (RFPN) and temporal noise (TN) as a function of well fill, rulingout the defective pixels detected by our "advanced" algorithm (black) or these defective pixels and all the RTS pixels. The fact that the curves are superimposed proves that the RTS pixels that have an influence on the RFPN are all detected by our "advanced" algorithm.

  defective + all RS pixels excluded

Table 2 -Criteria used to classify one pixel as a defective one (advanced algorithm).

 2

Proc. of SPIE Vol. 11002 110020E-2 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-3 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-4 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-5 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-6 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-7 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-8 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-9 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-10 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-11 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-12 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use

Proc. of SPIE Vol. 11002 110020E-13 Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 11 Oct 2019 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use


domain: spi phys

text:
P Catalano

D De Rosa

B Mele

R Tognaccini

F Moens

Effects of riblets on the performances of a regional aircraft configuration in NLF conditions

The application of the riblets to a typical regional turbo-prop configuration is discussed in this paper. The effect of the riblets is modelled as a singular roughness problem by a proper boundary condition proposed by University of Napoli. The model is described and first its application to 2D flows is presented. Then a typical wing-body of a regional aircraft is considered. The configuration has been designed to have a natural laminar flow for a large extent of the wing in cruise conditions. Riblets are applied at flow specifications representative of cruise in combination with the laminar flow technology and in climb/descent conditions. The save of fuel achievable in a typical operative day of such an aircraft has been also estimated.

I. Introduction

Today there is renewed interest in drag reduction mechanisms. Indeed, current concern over environmental pollution is forcing manufacturers to reduce pollutant emissions not only in the industrial field but also in the transport sector. Several methods for reducing viscous drag have been studied during the last three decades. These are used to delay laminar-boundary layer transition or to modify the turbulence structures of the boundary layer. Some methods are active controls, and others are passive as they do not require any energy use. Focusing on passive devices, the natural laminar flow technology aims at extending the laminar region as much as possible. Riblets consisting of stream-wise grooved surfaces are instead good candidates to reduce friction drag in the turbulent part of the flow.

A detailed and accurate simulation of the flow around riblets requires the adoption of Direct Numerical Simulations (DNS) or Large Eddy Simulations (LES). Neverthless, complex configurations at high Reynolds numbers can still only be studied by methods based the RANS equations. However, the RANS methods cannot resolve the scales of riblets (microns in aeronautical applications) and the effect of the riblets needs to be modelled. Mele and Tognaccini 1 have introduced the idea to model the effect of the riblets by a proper boundary condition at the wall such as usually done for the effect of roughness. The proposed boundary condition has been impelmented in two URANS flow solvers, the CIRA in-house developed UZEN code, and FLOWer, a code developed at DLR and used by University of Napoli. A relevant number of 2D and 3D test-cases have been carried out in order to validate and evaluate the model. [START_REF] Mele | Effect of riblets on a complex configuration in transonic conditions[END_REF][START_REF] Mele | A unified model for riblet simulation in complex flows[END_REF][START_REF] Mele | Analyses of non local effects of riblets[END_REF][START_REF] Mele | Performance assessment of a transonic wing-body configuration with riblets installed[END_REF] This paper deals with the effect of the riblets on the performances of a real aircraft when combined with the natural flow technology. Numerical simulations have been performed at different flow conditions representative of the characteristics points of a typical mission of a regional aircraft. The main goal is to evaluate the gain achievable by the riblets over an entire mission of the aircraft.

The paper is structured as follows. The modelling of the effect of the riblets is discussed in section II, and the results are presented in section III. First some results for 2D test cases are shown and then the application of riblets to the regional aircraft is discussed. The riblets are discussed in conjunction with the natural laminar flow technology. The weight that a new-generation turboprop aircraft could save is estimated. The conclusions are drawn in section IV.

II. Modelling the Effect of Riblets

It has been shown that riblets induce a displacement of the turbulence eddies with respect to the mean flow leading to a reduction of the momentum transfer to the wall with a consequent drag reduction. [START_REF] Bechert | Experiments on Drag-reducing Surfaces and their Optimization with an Adjustable Geometry[END_REF][START_REF] Luchini | Resistance of Grooved Surface to Parallel Flow and Cross-flow[END_REF][START_REF] Jiménez | Turbulent Flows over Rough Walls[END_REF] This mechanism is confined very near to the surface and induces a drag increase if riblets are large enough to interact directly with the turbulence structures. In practice, the effect of riblets is the shift of the origin of the velocity profile, i.e. a shift of the constant of the well known logarithmic law of wall:

U + = 1 κ a log(y + ) + B -∆U + , (1)
where the superscript + specifies non-dimensional quantities obtained by using wall variables, κ a is the Kármán constant, and B measures the influence of the wall on the velocity profile. Equation 1 is the same formula describing the effects of wall roughness on turbulent flows. The difference is that, in the case of roughness, ∆U + is usually positive returning an increase of drag, while riblets provide negative values of ∆U + with a reduction of the skin friction. Tani 9 re-analyzed the classical experimental data of Nikuradse 10 on turbulent flows over rough walls. He focused his attention on the transitional roughness regime that is characterized by a roughness height not large enough for a fully rough behaviour (usually the transitional roughness is defined for non-dimensional roughness height k + s < 50). He realized that sand roughness induces a reduction of skin friction if k + s , which is often defined as roughness Reynolds number, is lower than about 6. Above this critical value, the sand roughness induces an increase of drag. Tani already noted that this behaviour is surprisingly similar to riblet effects which however induce a much greater drag reduction. On the other hand, he emphasized that, unlike riblets, sand roughness is insensitive to flow direction.

The well known wall boundary condition of Wilcox for the k -ω turbulence models family [START_REF] Wilcox | Turbulence Modeling for CFD-II edition[END_REF] is widely adopted for simulating rough walls. Wilcox derived his boundary condition from Saffman [START_REF] Saffman | A Model for Inhomogeneous Turbulent Flow[END_REF] who first observed that, approaching the surface, the differential equations of the κ -ω turbulence models possess a kind of solution that returns a finite value for ω at the surface. This value can be, in principle, any value that correctly reproduces the logarithmic law of smooth, rough or grooved wall. Saffman proposed the following boundary condition for ω at the wall:

ω = ρu 2 τ µ • S R (k + s ) = τ w µ • S R (k + s ). ( 2
)
where ρ is the density, u τ = τ w /ρ is the friction velocity (τ w is the wall shear stress) and µ is the dynamic viscosity. Using equation 2 and from the analysis of the experimental results of Nikuradse, Wilcox derived his condition for S R on rough walls:

S R = (50/k + s ) 2 k + s < 25, 100/k + s k + s 25. ( 3
)
Thus ω at the wall depends on ∆U + or equivalently on S R = S R (∆U + ) and depends on the nature of the wall. In the case of transitional roughness S R (∆U + ) has been obtained verifying, by numerical experiments, the correct value of S R providing the ∆U + measured by Tani. A simple bell shaped function has been used to fit the new obtained values for S R :

S R = C 1 (k + s -C 2 ) 2n + C 3 , (4)
where C 1 = 1.666 • 10 3 ; C 2 = 3.25; C 3 = 0.1; n = 1. Equation 4 recovers for k + s = 0 and k + s = 6.5 the S R value for smooth wall.

The same equation 4 with different coefficients and written as a function of l + g = A + g , where A + g is the riblet non dimensional cross section area, can be used to simulate the effects of an arbitrary riblet family. In fact, following Mayoral and Jiménez, 13 l + g gives a better characterization of riblet performance than the non-dimensional riblet spacing s + and height h + (figure 1(a)). They found that the adoption of l + g collapse (a) Different shapes of riblets (b) "v-shape" riblets Fig. 1: Shapes of riblets riblet experimental data into a compact group of similar curves. Thus, in the case of riblets, equation 4 can be rewritten as:

S R = C 1 (l + g -C 2 ) 2n + C 3 . (5)
The coefficients have been obtained by numerical experiments matching the experimental data reported in: [START_REF] García-Mayoral | Drag reduction by riblets[END_REF] C 1 = 2.5 • 10 8 ; C 2 = 10.5; C 3 = 1. • 10 -3 ; n = 3. C 2 is equal to the value of l + g corresponding to the maximum value of S R , while C 1 and C 3 are related to the maximum value of S R . In the case of the riblet family reported by Walsh, 14 manufactured by 3M company, with a symmetric V-grooved section (figure 1(b)), the relation between l + g , s + and h + is s

+ = h + = √ 2 l + g .

III. Results and Discussion

A. 2D Flows

The model described in the previous section has been applied to a series of 2D flows for purposes of validation.

The results obtained are consistent with theoretical considerations and in good agreement with the available experimental data.

Flat Plate

A flat plate at Reynolds number (based on the length of the plate L) of 1 × 10 6 has been considered. The κ -ω SST turbulence model has been applied. The boundary condition 5 has been applied on all the length of the plate to model the effect of the riblets.

The drag reduction, computed as percentage variation between the drag coefficient of the riblets-on configuation C R D and the riblets-off configuration C S D , as function of l + g is shown in the left plot of figure 2. The maximum reduction is about 6% and is achieved, as expected, for l + g = 10.5. The same C D as the smooth (l + g = 0.) configuration is obtained for l + g between 20 and 22. Riblets with l + g = 22 provide an increase of drag.

The friction coefficients obtained by applying riblets with different l + g (according to equation 5) are shown in the right plot of figure 2. The Blasius laminar and Prandtl turbulent curves are reported as a reference. The riblets have been applied on all the surface of the plate also in the laminar part. No influence of the riblets is noted in the laminar zone of the flow but the numerical transition and the turbulent region is affected by the presence of the riblets. The transition is delayed by the riblets, and lower levels of the skin friction coefficient are attained in the turbulent part as l + g increases. The minimum level of C f is obtained, for l + g = 10.5 which is also the case that presents the most delayed transition. The friction coefficient is always lower with respect to the clean configuration except for the riblets with l + g = 22. According to the ∆C D results, the riblets provide the same C f as the clean configuration for 20 < l + g < 22. As discussed in the previous chapter (equation 1), the effect of the riblets can be seen as a shift in the log-law of the velocity profiles. This can be retrieved in the achieved results. Left plot of figure 3 shows the the theoretical curve (U + = y + ) is perfectly followed, while the logarithmic part of the profile is shifted by the riblets. The ∆U + due to the presence of the riblets can be estimated. In fact, at the edge of the boundary layer, the

∆U + E is ∆U + E = (U + E ) R -(U + E ) S = 1 κ a log(δ + ) R - 1 κ a log(δ + ) S + ∆U + (6)
where the subscripts S and R refer to the smooth and to the configuration equipped with riblets. By assuming that δ + is the same and considering that for the incompressible flow over a flat plate

U E ≈ U ∞ since 1 - U E U ∞ 2 = 2γ γ -1 P E -P ∞ ρU 2 ∞ (7)
the velocity shift is obtained as

∆U + = 2 C f R - 2 C f S ( 8
)
Right plot of figure 3 reports ∆U + for the turbulent region of the flat plate. As expected, it results to be greater than zero for all the values of l + g except for the case l + g = 22. The maximum achieved values of ∆U + are around 0.7 and is obtained for l + g = 10.5. The increment attainable for ∆U + is not linear with l + g but increases around the optimal value of 10.5. It is also interesting that this velocity shift is almost constant with the Reynolds number.

CAST7 Airfoil

The transonic flow around the CAST 7 airfoil is discussed in the following. Experimental data [START_REF] Coustols | Synthesis of Experimental Riblet Studies in Transonic Conditions[END_REF] The configuration has been considered with the transition imposed at 5% of the chord on both the upper and lower surfaces and the κ -ω SST turbulence model has been applied. The riblets are installed from the 15% of the chord.

The drag coefficient obtained for the clean configuration and with riblets of h = 0.023 mm. and h = 0.051 mm. is plotted as a function of the Mach number in figure 4. The results achieved are in a reasonable agreement with the experimental data. The behaviour with the Mach number is well predicted, but the C D is over-predicted for the smooth configuration (black color) and under-predicited for the riblets-on configurations with respect to the experiments (full symbols). In particular, the drag coefficient decreases at all the Mach numbers for the riblets with h = 0.023 mm. (red color) in both experimental and numerical data. Instead the riblets with h = 0.051 mm. present an increase of C D in the experiments, while the numerical simulations provide a slight decrease of drag at Mach 0.65 and 0.70 and a slight decrease at Mach 0.76. It turns out that the model and the numerics employed in the simulations require riblets somehow higher than in the experimets to achieve an increment of drag. The height of the riblets should be changed in the computations in order to get the same drag coefficient as the experiments. In particular the increase of drag presented in the experiments with riblets of h = 0.051 mm. is obtained in the numerical simulations with a height of 0.064 mm. at Mach 0.76, h = 0.059 mm. at Mach 0.70, and h = 0.056 at Mach 0.65.

As an example, the friction coefficient for the smooth configuration and for the configuration with riblets of 0.023 mm. height is presented in left plot of figure 5. The decrease of the friction coefficient is evident The physical height of a symmetric "v-shape" riblet can be related through to the viscous height h + through the following relation:

h + = hu τ ρ µ = h L REF Re ∞ C F 2 ρ ρ ∞ µ ∞ µ (9)
The above equation 9 is useful to evaluate the physical height of "v-shape" riblets corresponding to the "optimal" h + . This is shown in right plot of figure 5 at the three Mach numbers taken into consideration.

The height increases in a not linear way with the Mach number. It also increases along the chord on both the upper and lower surface of the airfoil and presents a spike, almost doubling, close to the shock wave.

DU 96-W-180 airfoil

Chamorro et al. [START_REF] Chamorro | Drag reduction of large wind turbine blades through riblets: evaluation of riblet geometry and application strategies[END_REF] tested various riblet size and shapes provided by 3M company on a DU 96-W-180 airfoil.

The experiments were performed in low speed conditions at a Reynolds number (based on the chord of the airfoil) of 2.2 × 10 6 and varying the lift coefficients from 0.25 to 1.14.

The same airfoil was studied in the experiments of Sareen et al. [START_REF] Sareen | Drag reduction using riblet film applied to airfoils for wind turbines[END_REF] They also tested various riblet size provided by 3M company at three different Reynolds numbers 1.0 × 10 6 , 1.5 × 10 6 and 1.85 × 10 6 .

The effect of full or partial riblet coverage on the airfoil has also been investigated in both experiments. The authors calculated, by Xfoil simulations, the transition to turbulence over the clean airfoil in its operational range and then, in the case of partial riblet coverage, the riblet film was applied starting from 40% of the airfoil chord on the upper surface and 70% on the lower surface. Sareen et al. 17 also determined the experimental separation and reattachment points and reported an oil flow visualization of the detected laminar separation bubble at α = 6 • and Reynolds number 1.5 × 10 6 .

The present numerical simulations have been performed at the same flow conditions as the experiments with partial riblet coverage. All calculations have been carried out at M ∞ = 0.1 and specifying the transition on the upper and lower surface at the same position reported in the experiments. In the laminar zone the production terms of the turbulence equations are set to zero; such transition treatment guarantees optimal flow prediction in case of laminar separation bubbles. [START_REF] Catalano | On the implementation of a turbulence model for low Reynolds number flows[END_REF] The computed drag reduction vs. riblet height at Re ∞ = 1.5 × 10 6 is compared with the measurements of Sareen et al. [START_REF] Sareen | Drag reduction using riblet film applied to airfoils for wind turbines[END_REF] The computed skin friction coefficients with and without riblets at Re ∞ = 1.5×10 6 and C l = 0.75, 0.875, 1.0 are shown in fig. 8. Detected laminar separation bubbles are clearly visible and the separation and reattachment points are in good agreement with experimental data. It is worth noting that the whole effect of riblets on skin friction coefficients is in the attached flow zone. The skin friction seems to be not influenced inside the laminar separation bubble .

B. Regional Aircraft Configuration

A wing-body configuration of a typical regional aircraft (figure 9) has been taken into consideration. [START_REF] Catalano | Influence of riblets and laminar flow control on the performances of a regional aircraft[END_REF] The wing has been designed following the process used for a transonic application. [START_REF] Salah | Natural laminar flow transonic wing design applied to future innovative green reegional aircraft[END_REF] The airfoil sections have a laminar flow in cruise conditions and the twist is optimized for the low-speed performances. A structured Different assumptions are made in order to assess the effect of the riblets in combination with a natural laminar flow (NLF). Riblets are applied with the assumption of "fully turbulent" flow and in case the flow is assumed laminar for a large part of the wing according to the NLF technology.

The riblets have been installed only on wing, only on body and on all the configuration. This kind of analysis allows to determine the zones of the configuration where the riblets are more effective.

The flow conditions investigated are representative of the mission typical of a regional aircraft. The main aim is to evaluate the gain in terms of aerodynamic drag that could be achieved by applying the riblets to such a configuration.

Cruise conditions with fully turbulent assumption

A flow condition corresponding to Mach number 0.50 and Reynolds number (based on the chord of the wing) of 3.0 × 10 7 has been considered.

The ∆C D (pressure and friction contributions) with respect to the clean configuration achieved by applying riblets with l + g = 10.5 is reported in figure 10. Riblets only on wing are more effective at the low incidences, while riblets only on body provide a gain in drag quite constant with α. A combined effect is obtained when riblets are applied over all the configuration.

In percentage terms, the maximum gain in total drag evaluated as:

∆C D = C R D -C S D C S D
is about 9% and decreases with the lift coefficient. The gain in friction drag evaluated as

∆C Df = C R Df -C S Df C S Df is constant with
the lift coefficient and is about 11%. Figure 11 reports the skin friction distribution on the upper surface of the configuration in riblets-on and riblets-ff conditions. The riblets have been installed only on the wing. A clear decrease of the skin friction is obtained in the front part of the wing.

The right plot of figure 11 reports, as an example, the C f at a station located at 40% of the wing span. The gain in C f decreasing as x increases is clearly visible. The riblets-on and clean configuration have the same levels of friction in the trailing-edge region.

The effect of applying the riblets only on the body can be appreciated in figure 12 The zone of the body where the riblets are more effective are the nose and wing-fuselage intersection.

The right plot of figure 12 shows the friction distribution along a constant-y section of the body. An effect of riblets, not visible in figure 12(a), is also present in the rear part of the fuselage.

Cruise conditions with natural laminar flow technology

The wing of the configuration has been designed to have a laminar flow for a large extent of the wing in cruise conditions. This assumption has been numerically reproduced by imposing a transition at 50% of the local chord in the outboard zone of the wing. The inboard region has been assumed turbulent because the presence of the engines prevents to keep the flow laminar (figure 13). The effect of applying the natural flow technology can be appreciated in figure 14 where the comparison between the "fully turbulent" (with and without riblets) and NLF configurations is shown. The ∆C D is evaluated with respect to the clean "fully turbulent" configuration. It is possible to note how the order of magnitude of ∆C D achieved by the NLF technology is the same as for riblets applied on both wing and body on "fully turbulent" configuration (figure 14a). In percentage, the maximum ∆C N LF D is about 12% and decreases with the α up to about 2.5%. The gain for the friction component of the drag is about 12% and keeps constant with the incidence.

The NLF condition presents a large region of turbulent flow and therefore the riblets can be effective also for such a configuration. The right plot of figure 14 shows the ∆C D achieved applying the riblets together with the NLF technology. In this plot, the ∆C D is evaluated with respect to the clean configuration with the NLF assumption. In case of riblets on both body and wing, a further decrease of the total drag coefficient varying with α from about 9% (3% body and 6% wing) to about 3% (2.5% body and 0.5% wing) has been achieved. The friction drag decreases of about 12% (8% body and 4% wing).

Finally the contour plot of the friction coefficient over the upper surface of the NLF configuration in riblets-off condition and with the riblets applied on both body and wing is presented in figure 15. The effect of the riblets is visible on the nose of the body and in the zone of the intersection between the fuselage and the wing. The C f at an outboard station of the wing is presented in the right plot of the figure 15. It is possible to note how the NLF technology presents in the laminar zone with respect to the "fully turbulent" assumption a reduction in friction larger than the riblets. However, the riblets keep their effectiveness in the turbulent region and provide a further gain in C f with respect to the NLF technology.

Climb/Descent condition

A flow specification of Mach number 0.20 and Reynolds number (based on the chord of the wing) of 1.30×10 7 has been chosen as representative of a climb/descent condition. Only "fully turbulent" conditions have been taken into consideration.

The gain in drag coefficient in climb conditions is shown in figure 16. The "optimum" riblets have been applied on wing, on body, and on all the configuration. Riblets on wing have a similar effect, in term of ∆C D as the riblets on body. The maximum gain in drag is achieved when riblets are applied on wing and fuselage and results to be about 7% and decreases with the incidence. The ∆C Df is about 8% and is constant with the α.

Finally, the contour plot of the friction coefficient over the upper surface is presented in figure 17. Again, the riblets are effective in the front part of the wing with an effect quite constant in the span-wise direction and in the region of the nose of the body and of the wing-body intersection.

Estimate of gain for an operative day

The Breguet formula for the range D . This is obtained by considering the drag breakdown of a typical turboprop-driven aircraft as shown in figure 18b where "other" is represented by interferences, windshield,

R = η gc s C L C D ln W i W f (10)
C D = C W B D + C HT D + C V T D + C N a D + C Ot D = C W B D 1 + 0.059 + 0.073 + 0.176 + 0.162 (13)
The drag coefficient, evaluated by equation 13 by exploiting the drag polars computed for the wingbody configuration, is used in equation 10 where η is the efficiency of the engine and c s is the specific fuel consumption. The values considered typical for a turbo-propeller driven aircraft are η = 0.85 and c s = 0.70(lb/hph). The difference in weight ∆W that can be achieved by riblets and NLF technology has been computed with respect to W c f , the weight at the end of the cruise phase of the clean configuration. The ∆W , in percentage terms, that a turboprop regional aircraft could save in a design mission with a cruise of 900 Nm. has been estimated and reported in table 1. The weight that a turboprop regional aircraft The analysis above reported is partial and suffers from some assumptions especially with regards to the evaluation of the drag coefficient. However, it is possible to state that the NLF technology would allow to save more weight than riblets, about 30% more. Is is also worth noting that, in terms of saving weight, riblets on wing are more effective (almost two times) than riblets on fuselage.

IV. Conclusions

The application of the riblets, a promising device for the drag reduction, to a new-generation turboprop configuration has been presented. RANS simulations have been performed by exploiting a proper boundary condition to simulate the effect of the riblets. Drag polars in cruise and climb/descent conditions have been computed.

The wing of the configuration has been designed with the aim to have laminar flow in cruise conditions, and the riblets have been considered also in conjunction with the NLF technology.

The weight that a turbopropeller-driven aircraft could save in the cruise phase of typical mission by applying these drag-reduction devices has been estimated.

NLF technology is shown to be more effective but performance improvement by riblets is comparable and significant.

V. Ackowledgements

The research leading to these results has received funding from Clean Sky 2 Joint Undertaking under the European Union's Horizon 2020 research and innovation programme under grant agreement CS2-REG-GAM-2014-2015-01.

Fig. 2 :

 2 Fig. 2: Flat plate at Re = 1 × 10 6 : Drag reduction and Friction coefficient

Fig. 3 :

 3 Fig. 3: Flat plate at Re = 1 × 10 6 : Effect of riblets on velocity

  Fig. 4: CAST-7 Airfoil at Mach = 0.76, Re = 3 × 10 6 , and α = 0:Drag Coefficient. Full Symbols: Experimental; Void Symbols: Numerical.

Fig. 5 :

 5 Fig. 5: CAST-7 Airfoil at Re = 3 × 10 6 , and α = 0: Skin friction coefficient and height of the riblets

  in fig.6. A satisfactory agreement is shown. The achieved maximum drag reduction is (a) Re∞ = 1.5 × 10 6 . : C l = 0.75, •: C l = 0.875, : C l = 1. (b) C l = 0.75. •: Re∞ = 1.5 × 10 6 , : Re∞ = 1.85 × 10 6 , : Re∞ = 1 × 10 6 .

Fig. 6 :

 6 Fig. 6: DU 96-W-180 airfoil. Computed and experimental drag reduction vs. riblet height. Void symbols: Numerical, Full Symbols: experimental

Fig. 7 :

 7 Fig. 7: DU 96-W-180 airfoil with symmetric V-grooved riblets of h = s = 80µm, Re ∞ = 2.2×10 6 . Computed and experimental drag reduction vs. lift coefficient. -∆-: computed, : experiments 16

  (a) C l = 0.75. (b) C l = 0.875. (c) C l = 1.0.

Fig. 8 :

 8 Fig. 8: DU 96-W-180 airfoil, Re ∞ = 1.5 × 10 6 . Computed skin friction coefficients. --: computed with riblets h = 62µm, ---: computed without riblets, : experimental separation and reattachment points without riblets.

Fig. 9 :

 9 Fig. 9: Regional aircraft configuration

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Wing-body configuration: ∆C D (percentage gain) by riblets. --: ∆C D , ---:∆C Df

  Fig. 12: Wing-body configuration: Skin friction coefficient for riblets only on body

Fig. 13 :

 13 Fig. 13: Wing-body configuration: Numerically reproducing the natural laminar flow technology. blue region: laminar; red region: turbulent

Fig. 14 :

 14 Fig. 14: Wing-body configuration: NLF technology and riblets. --: ∆C D , ---:∆C Df

Fig. 16 :Fig. 18 :

 1618 Fig. 15: Wing-body configuration: Skin friction with NLF technology and riblets

Table 1 :

 1 Saving of weight for a design mission with cruise phase of 900 Nm.Clean Riblets wing Riblets body Riblets wing and body NLF NLF+Riblets wb 100 * ∆W/W c a typical mission with a cruise of 350 Nm. has been estimated and reported in table 2.

	f	0	3.75	1.74	7.38	9.86	16.03
	could save in

Table 2 :

 2 Saving of weight for a typical mission with cruise phase of 350 Nm. Clean Riblets wing Riblets body Riblets wing and body NLF NLF+Riblets wb 100 * ∆W/W c

	f	0	1.6	0.74	3.15	4.20	6.82