Cédric Huchette
email: cedric.huchette@onera.fr

Gillian Leplat
email: gillian.leplat@onera.fr

Valentin Biasi
email: valentin.biasi@onera.fr

DAMAGE ANALYSIS FOR THERMAL LOADING INDUCED BY LASER IMPACT IN EPOXY COMPOSITE LAMINATE

Keywords: fire, damage, delamination, thermal, modelling

The thermal response of a composite material subjected to high heat flux is driven by heat conduction combined with endothermic (pyrolysis) and exothermic (oxidation) chemical reactions that change locally the thermo-physical properties and generate decomposition gas flowing throughout the porous charring medium. Heat and mass transfer is also affected by mechanical damage experienced by laminates such as matrix or interfacial cracking the origin of which may be related to internal stresses and thermal loads as well as pore and char formation. The challenge is to analyse concurrently the multi-physics phenomena in order to identify and understand the damage processes.

Introduction

The thermal response of a composite material subjected to high heat flux is driven by heat conduction combined with endothermic (pyrolysis) and exothermic (oxidation) chemical reactions that change locally the thermo-physical properties and generate decomposition gas flowing throughout the porous charring medium. Heat and mass transfer is also affected by mechanical damage experienced by laminates such as matrix or interfacial cracking the origin of which may be related to internal stresses and thermal loads as well as pore and char formation. The challenge is to analyse concurrently the multi-physics phenomena in order to identify and understand the damage processes.

However some progress has been made in understanding the composite thermo-mechanical behaviour when subjected to a fire [START_REF] Mouritz | Fire properties of polymer composite materials[END_REF][START_REF] Mouritz | Review of fire structural modelling of polymer composites[END_REF], it remains the combination of complex interactions between the condensed phase and the gas phase within the material but also at the interface with the flame. It is all the more important since fire dynamics considerations are mandatory to assess relevant space and time distribution of the heat exchanges from the flame to the material surface. The conventional configurations and experimental setups [START_REF] Mouritz | Fire properties of polymer composite materials[END_REF][START_REF] Mouritz | Review of fire structural modelling of polymer composites[END_REF] are designed to limit the analysis to the through-thickness direction while major composite laminates have an intrinsic orthotropic behaviour caused by the fibre reinforcement and the stacking sequence for laminates. The in-plane behaviour is unfortunately avoided while it may prevail in a practical event.

The present study is devoted to the analysis of the thermo-chemical phenomena and the induced damage occurring within decomposing aeronautical composite materials. A dedicated experiment is thus proposed where accurate test conditions are achieved, inhomogeneous heating is applied, nonintrusive transient measurements are performed and fire dynamics discrepancies are prevented. The material thermal response is correlated with micrographic imaging carried out post-decomposition to investigate how and why delamination damage occurs and grows.

Materials and methods

Material

The investigation deals with a material used in the aeronautical industry for primary and secondary aircraft structures. The T700GC/M21 is a composite laminate consisting of TORAY T700GC carbon fibres and HEXCEL M21 epoxy resin reinforced by thermoplastic nodules. Sixteen 260 µm -thick M21/35%/268/T700GC unidirectional 268 𝑔. 𝑚 ! prepreg plies are stacked and cured to manufacture the composite laminate. This high fibre mass fraction is the consequence of the stacking of two fibre tows in each prepreg tape. The M21 resin is a tough epoxy matrix and exhibits excellent damage tolerance, especially in the case of high energy impacts. To obtain such properties, it contains the thermoplastic nodules around each fibre tow with a characteristic length of 25 µm (Fig. 1). The fibres average diameter is 7 µm and the volume fraction of fibres is 0.57 for the cured material. A 16 ply quasi-isotropic laminate with a final material thickness of 4.16 mm has been studied with the following stacking sequence [0/45/90/-45/0/45/90/-45] s .

Figure 1. Micrography of the virgin material. T700GC/M21 intra and interply interfaces clearly visible due to the presence of thermoplastic nodules and to the stacking of two fibre tows within each ply.

BLADE facility

The BLADE facility (characterisation and decomposition laser facility from the French ''Banc Laser de cAractérisation et de DEgradation'') is devoted firstly to the characterisation of thermo-physical properties of anisotropic materials and secondly to the analysis of the thermal response when they decompose. Even if the apparatus is original, its principle is very simple (Fig. 2). The apparatus offers two main advantages. First, the enclosure is air-filled, temperature-regulated and the pressure within is decreased down to 3 mBar which is low enough to avoid any convective heat transfer. This last feature makes the BLADE facility original and it results in accurate control of test conditions. Indeed, the determination of the convective heat transfer coefficient would have been very difficult even in a simple geometrical configuration without any guarantee of a good representativity of the convective heat transfer within the chamber. The second advantage of the BLADE facility consists in using a laser instead of a flame to heat up a square-shaped, 80 mm-wide, test coupon during the experiment. The heat flux magnitude and distribution on the material exposed surface may be very difficult to assess accurately and reliably when using standard flames generated by a burner. Instead, the laser provides a heat source, stable in space and time, and it is characterized experimentally preliminary to the tests using inverse heat conduction method on a pure titanium reference material. The Gaussian distribution of the beam collimated at ∅ = 21.8 𝑚𝑚 at 1/𝑒 ! creates a non-uniform heating on the exposed surface which helps revealing the anisotropic behaviour of such composite laminates and offers ultimately valuable data for multi-dimensional models validation. The transient temperature at the back (unheated) side of the test coupon is measured using an infrared thermography camera from the test coupon at rest, then during the heating phase up to the cooling phase when the shutter closes its diaphragm and the laser is switched off.

Results and damage analysis

Thermal behaviour of decomposing laminates

The IR thermography measurements performed on the T700GC/M21 QI test coupons subjected to identical conditions are presented in Fig. 2. The instantaneous temperature maps extracted for each test at 𝑡 = 150 𝑠 are presented on the left hand side of Fig. 3. Temperature evolutions as a function of time are also extracted along a horizontal mid-height profile at different locations from the centre and plotted on the right hand side Fig. 3.

A gradual temperature rise of more than 250 K from room temperature is measured on the back surface. The thermal response remains axisymmetric during this first phase and the highest temperature magnitude is observed at the centre of the coupons where the thermal loading at the front side is the highest. Repeatability is evaluated from these 3 tests performed on 3 different coupons and all of them show similar thermal behaviour in the first part of the experiment up to t = 150 s (Figs.

34).

No in-plane deformation is detected from the IR measurements considering the camera spatial resolution.

The behaviour of the unexposed surface resulting from axisymmetrical Gaussian laser heating on the front surface can thus be addressed under a thermal loading high enough to cause chemical decomposition. The repeatable temperature rise observed in the first part of the experiment is governed by the energy balance between the incoming laser flux, heat conduction within the condensed phase and radiative exchanges with the temperature-regulated chamber walls. Then, temperature-activated chemical reactions occur and affect the resin through a pyrolysis charring process and probably also affect the new-formed char phase through an oxidation process due to the high temperature magnitude expected on the front face of the material. The temperature maps up to t = 150 s (Fig. 3) show that the thermal response of the QI laminates to an axisymmetrical centred thermal loading remains isotropic in the composite panel plane. It is important to note than no ignition is observed all along the experimental process.

After t = 150 s, the three coupons present for different time, sudden drops of temperature. These temperature drops are not observed for from the centre of the plate. These drops are the consequence of cracks apparition in the thickness of the laminate and there influence on the thermal response on the opposite face could be linked with the size of these defaults (less than 50 mm wide).

Post mortem analysis

In order to validate the hypotheses made for the thermal analysis, post-decomposition micrographic analysis is carried out after each test. The micrographs are obtained from longitudinal cross sections of the coupons. A diamond saw device is used and the observed edges are polished from coarse to fine grain polishing paper in order to obtain a mirror finish. Olympus microscope with three axes controller and a monochromatic digital camera (Sensicam PCO) are used to capture the 205 pictures required to observe a whole cross section. The total image reconstruction of the polished edges presented hereafter is performed with an in-house code. Longitudinal median cross-sectional views are presented for each coupon in fig. 5 in order to observe the post-test material condition. Similar observations are taken from transversal median cross sections and are not presented in this paper. For each coupon, the position of the charred zone (in black) and the major delamination cracks (red lines) are reported in a sketch representing the stacking sequence and the damage patterns. The charred zone is idealised by the shape of half an ellipse and the delamination cracks by straight lines. The absolute position of the crack tips for each major delamination crack is reported in Table 1. It is essential to note that only major delamination cracks are reported, although several matrix cracks are present around the charred zone and most of them are under this zone.

Thanks to these micrographs, several comments could be made. First of all, all coupons present delamination cracks, however, the number, theirs positions and their size are different for the 3 coupons. The number and the size of delamination cracks are correlated with the number and the depth of the observed temperature drops on the thermal response of the opposite face to the thermal loading (Fig. 4). For example, the biggest temperature drop is observed for the coupon #3 and its correlated with the biggest delamination crack noticed between the first plies of the laminate #3 (close to the opposite face to the thermal loading). On the opposite, the coupon #2, which presents several temperature drops, have three major matrix delamination cracks. Concerning, the charred zone, the depth and the wide of this zone are quite similar, nevertheless, we can observed for the coupon #2 that its depth one is bigger than the two others. The raison of this difference seems to be due to the apparition of a wide delamination crack in the middle of the laminate thickness. The delamination crack reduces the heat transfer in the through-thickness direction and acts as a thermal barrier, which protects the back surface from thermal decomposition. As a consequence, the earlier the delamination crack occurs, the thicker the charred area is because the front part keeps on decomposing while the back part remains temporarily isolated. Table 1. Location and length of the main delamination cracks for coupons #1-3 measured from micrographic images presented in Fig. 5 Coupon

Y l [m] Y r [m] Z l [m] Z l [m] L [m] 1 -13.03 ×10 -3 7.19 ×10 -3 1.85 ×10 -3 2.22 ×10 -3 20.22 ×10 -3 1 -6.95 ×10 -3 6.28 ×10 -3 2.54 ×10 -3 2.60 ×10 -3 13.23 ×10 -3 2 -17.51 ×10 -3 20.99 ×10 -3 2.49 ×10 -3 2.56 ×10 -3 38.49 ×10 -3 2
-11.16 ×10 -3 5.75 ×10 -3 3.75 ×10 -3 3.89 ×10 -3 16.91 ×10 -3 2 -12.43 ×10 -3 6.00 ×10 -3 1.78 ×10 -3 1.96 ×10 -3 18.42 ×10 -3 3 -21.42 ×10 -3 22.39 ×10 -3 3.34 ×10 -3 3.43 ×10 -3 43.81 ×10 -3 3 -6.37 ×10 -3 5.49 ×10 -3 2.01 ×10 -3 2.35 ×10 -3 11.87 ×10 -3

The delamination cracks experienced by composite laminates can be associated with both thermal and mechanical loadings. Thermally induced damage concerns the chemical decomposition of matrix and fibres, the char formation associated with gas production as well. Those phenomena are mainly responsible for the onset of diffuse delamination observed in the experiments. Indeed, the gas pores, resulting from the thermo-chemical decomposition of the matrix that occurs first as the temperature increases, facilitate the onset of microscopic cracks within the condensed phase. This material presents very large resin-rich regions at ply interfaces but also within each ply as illustrated in 1). The first delamination crack (in the red box) is located near the charred zone, at the intraply interface of the internal 90° ply. The matrix cracks are highlighted in yellow and can be observed on the right of the delamination crack tip, with very small matrix damage in the intraply interface of the same ply. A higher magnification examination of these small matrix cracks reveals charred thermoplastic nodules (in black in the yellow box). These observations support the assumption of thermally-induced delamination because the matrix reacts first to a temperature rise to create a solid residue (char) and gas pores acting as delamination precursors. Only the matrix (resin and thermoplastic nodules) pyrolysis is characterised by a noticeable mass loss in the thermo-gravimetric analysis performed on this material. These micrographs show that matrix decomposition starts with nodule decomposition. Moreover, they reveal that matrix decomposition can start far from the charred zone, due to elevated temperature. for estimate the visible char formation and imposes that the mass fraction of the polymer matrix is reduced by ∼ 20%, due to decomposition and vaporisation. In the Biasi's model [START_REF] Biasi | Modélisation Thermique de la Dégradation d'un Matériau Composite[END_REF], this state corresponds to a value of 0.2 for the pyrolysis reaction growth (α Pyro ). Moreover, these numerical results validate the assumption of thermal degradation "far" from the charred zone and the coupling between thermal degradation and mechanical loading for the onset of delamination cracks in the laminate. The implementation of the Biasi's model in a finite element code will permit in future work to estimate, by thermo-mechanical analysis, the stresses at the interface, to simulate the damage evolution and the thermal response of the material and to predict the residual strength of the laminate with the objective to design structural part.

Conclusions

The current study is based upon the analysis of the composite thermal response subjected to laserdriven heating within a test chamber. As a direct consequence, both experimental surrounding conditions and the heat source supplying the thermal loading distribution into the material exposed surface are accurately controlled and provide a relevant framework for understanding the material behaviour and for developing and validating representative models. The analysis of the transient thermal response with the support of the unsteady numerical simulation and valuable micrographic imaging performed post-decomposition has revealed that heat transfer from the exposed face makes the temperature of the condensed phase increase highly enough to cause pyrolysis of the matrix into a significant depth of the laminate. The orthotropic heat transfer induced by the fibre orientation and the stacking sequence, associated with the inhomogeneous distribution of the heat flux has helped in identifying the physical mechanisms linking thermal decomposition to delamination cracks onset and growth. Future work will focus on the estimation of the internal stress due to the material degradation and the thermal loading in order to estimate the damage scenario. Moreover, a coupling approach between the thermal solver and a thermo mechanical one will permit to estimate the influence of the matrix cracking on

Figure 2 .

 2 Figure 2. Principle of the BLADE facility.

Figure 3 .

 3 Figure 3. Thermal response on the unexposed surface measured by infrared thermography. Left: instantaneous temperature maps at t = 150 s. Right: temperature extracted over time at different locations from the centre of the coupon.

Figure 4 .

 4 Figure 4. Comparison between experimental results and numerical ones for the evolution of the temperature on the centre of the face opposite to the laser impacted face.

  (a) Sketch of the location and length of the main delamination cracks reported in Table1

3 Figure 5 .

 35 Figure 5. Micrographs of the 3 coupons longitudinal cross sections associated with sketches of the major delamination cracks and charred zone for the different specimens.

Fig 1

 1 due to the presence of these thermoplastic nodules. These resin-rich regions offer preferential paths for the cracks initiation and propagation. Fig 6 focuses on two examples of significant delamination cracks for coupon #2 (second and third major delamination cracks of coupon #2 in Table

Figure 6 .

 6 Figure 6. Influence of the charred nodules on the onset of a delamination crack (longitudinal cross section of coupon #2).

Figure 7 .

 7 Figure 7. Comparison between experimental (coupon #1 on left) and numerical estimation (on the right) of the charred zone.

Acknowledgments

This work is funded by ONERA within the frame of an internal collaborative project named MADMAX (French acronym for "Modèles Aérothermomécaniques de Dégradation de MAtériaux compleXes") dedicated to the development of aerothermomechanical models for decomposing complex materials.


Christopher R Yost

P Minnis

K M Bedka

L Nguyen

R Palikonda

D A Spangenberg

J W Strapp

J Delanoë

A Protat

Radar Measurements of IWC

Introduction

A23A-0185

At least one hundred jet engine power loss events since the 1990s have been attributed to the phenomenon known as ice crystal icing (ICI). Ingestion of high concentrations of ice particles into aircraft engines is thought to cause these events, but it is clear that the use of current on-board weather radar systems alone is insufficient for detecting conditions that might cause ICI. Passive radiometers in geostationary orbit are valuable for monitoring systems that produce high ice water content (HIWC) and will play an important role in nowcasting, but are incapable of making vertically resolved measurements of ice particle concentration, i.e., ice water content (IWC). Combined radar, lidar, and in-situ measurements are essential for developing a skilled satellite-based HIWC nowcasting technique. The High Altitude Ice Crystals -High Ice Water Content (HAIC-HIWC) field campaigns in Darwin, Australia; Cayenne, French Guiana; and Ft. Lauderdale, FL, have produced a valuable dataset of in-situ total water content (TWC) measurements with which to study conditions that produce HIWC. The NASA Langley Satellite ClOud and Radiative Property retrieval System (SatCORPS) was used to derive cloud physical and optical properties such as cloud top height, temperature, optical depth, and ice water path from multi-spectral satellite imagery acquired throughout the HAIC-HIWC campaigns.

These cloud properties were collocated with the in-situ TWC measurements in order to characterize cloud properties in the vicinity of HIWC. Additionally, a database of satellite-derived overshooting cloud top (OT) detections was used to identify TWC measurements in close proximity to convective cores likely producing large concentrations of ice crystals. Certain cloud properties show some sensitivity to increasing TWC and a multivariate probabilistic indicator of HIWC was developed from these datasets. This paper describes the algorithm development and demonstrates the HIWC indicator with imagery from the HAIC-HIWC campaigns. Vertically resolved IWC retrievals from active sensors such as the Cloud Profiling Radar (CPR) on CloudSat and the Doppler Radar System Airborne (RASTA) provide IWC profiles with which to validate and potentially enhance the satellite-based HIWC indicator. All other parameters except effective ice particle radius R eff show some ability to distinguish LIWC from HIWC as indicated by the separation between the green and black curves. However, it is difficult to distinguish HIWC from moderate values (small separation between the blue and green curves).

In-situ Measurements of IWC

Overview of HAIC-HIWC Campaigns and Datasets

Effective ice particle radius R eff exhibits no skill for distinguishing even very low from very high IWC.

IWP is a function of both COD and R eff , so it is preferable to use COD rather than IWP as a HIWC predictor.

TWC is a flight-level measurement.

RASTA radar retrievals of IWC typically reveal higher IWC below flight level.

Spatial/temporal collocation of aircraft and satellite data:

Assuming a typical airspeed of 180 m/s, the 1-Hz aircraft measurements were averaged to 45-s intervals in order to approximate the relatively coarse satellite spatial resolution. The four nearest satellite pixels were matched to each segment of the flight tracks and the mean cloud properties were computed. Temporal differences between the observations were restricted to less than 10 minutes. Satellite temporal resolution ranged from 1-minute (GOES-14 imagery for three flights during Florida campaign) to 30-minute (GOES-13 during Cayenne campaign).

Summary

High concentrations of ice particles, i.e., HIWC, have been linked to aircraft engine power loss events at cruise altitudes where the existence of supercooled water droplets is extremely improbable.

In-situ ice water content measurements were matched to satellite cloud property retrievals in order to characterize HIWC events from a satellite perspective and develop an algorithm to estimate the probability of high ice water content PHIWC in satellite imagery.

Most HIWC occurred within or near optically thick clouds with significant texture and/or prominent cold spots. There is significant overlap in the distributions of different satellite observations for moderate and high IWC encounters.

PHIWC validation studies highlight the tradeoff between probability of detection and false alarm rates. The HIWC nowcasting community will need to establish a method to effectively validate HIWC diagnostic products.

Next-generation sensors such as the Advanced Baseline Imager (ABI) on GOES-R and the Advanced Himawari Imager (AHI) on Himawari-8 offer the necessary spatial and temporal resolution for studying environments that produce HIWC.

  Images (derived from MTSAT-1R) of (a) visible reflectance, (b) IR window brightness temperature (BTW), (c) visible texture rating (magenta contours) and OT probability rating (color shading), (d) difference between IR-window and tropopause temperatures (ΔBTW), (e) brightness temperature difference between water vapor and IR-window channels (BTD), (f), distance to nearest OT, (g) ice cloud optical depth (COD), and (h) ice particle effective radius for Darwin flight #22 on 18 February, 2014." During the HAIC-HIWC campaigns, in-situ cloud total water content (TWC) was measured by an isokinetic evaporator probe (IKP2; Davison et al. 2012, AIAA) mounted to a research aircraft. Satellites in geostationary orbit were relied upon for continuous multi-spectral imagery of storms producing large concentrations of ice particles. A brief summary of the HAIC-HIWC campaigns is provided in the table below including the geostationary satellite imagers that provided coverage for each campaign. Cloud property retrievals (e.g., cloud boundaries, optical depth, water path) were derived from satellite imagery using the SatCORPS (Minnis et al. 2011), and an automated overshooting top (OT) detection algorithm (Bedka and Khlopenkov, 2016, JAMC) was used to identify locations of cirrus anvils and deep convective updrafts.The OT detection algorithm is a set of statistical, spatial and frequency analyses of IR and VIS imagery and consists of two main components:(1) IR-based detection of anvil clouds and brightness temperature (BTW) minima within these regions.OT candidates are assigned an "OT Probability" score.(2) Quantification of cloud texture in 1-km visible imagery via a unitless rating. Values greater than 5 are indicative of strong vertical motions and gravity waves, and larger values correspond to classic OT "cauliflower-like" texture.

  No single threshold can delineate HIWC from low IWC (LIWC), but values smaller than 0.25 g m -3 are typically not considered high. Values of 0.25 and 1.0 g m -3 represent the 48 th and 82 nd percentiles, respectively, of the TWC* dataset (45-s averages). Cumulative distribution functions (CDFs) of satellite observations/retrievals were generated for different ranges of TWC (below). Distance from OT (d OT ) shows the most potential for distinguishing LIWC from HIWC. Almost all instances of TWC > 1.0 g m -3 occurred within ~40 km of evident convection identified by the automated OT detection algorithm. **Delanoë, J., 2016. Radar Aéroporté et Sol de Télédétection des Propriétés Nuageuses. (RASTA), Version 5.2. Accessed 28 July, 2016. *Strapp, W. 2014. French Falcon IsokineQc Evaporator Probe (IKP2) Data, Darwin, Version 5.0. UCAR/NCAR --Earth Observing Laboratory. hZp://data.eol.ucar.edu/dataset/384.011. Accessed 4 Nov 2016. *Strapp, W. 2016. French Falcon IsokineQc Evaporator Probe (IKP2) Data, Cayenne, Version 5.0. UCAR/NCAR --Earth Observing Laboratory. hZp://data.eol.ucar.edu/dataset/486.001. Accessed 22 Aug 2016. CDFs of cloud properties for three ranges of IKP2 TWC. Figure includes Florida IKP2 draft dataset.! Mean TWC was computed as a function of ΔBTW, COD, and d OT (left). Since visible texture detection is not available at night, mean TWC was computed as a function of distance from IR-only detections too. A fuzzy logic method was adopted here for estimating HIWC probability PHIWC) in satellite imagery given ΔBTW, COD, and d OT . Functional fits to mean TWC (left, magenta curves) were used to obtain values ranging 0.0 -1.0 for each input. Final PHIWC values were obtained by taking the cubic root of the product of the three scores. TWC and PHIWC trends track especially well with trends in d OT , particularly when the satellite temporal resolution is very high (see time series below). During the Florida campaign, 1-minute imagery was available for three flights (see right time series). RASTA** (Doppler Radar System Airborne) is an airborne 95-GHz cloud radar operated on the Falcon-20 during the Darwin and Cayenne campaigns. RASTA radar measurements are capable of providing IWC retrievals for the entire cloud column thus providing a more complete characterization of cloud structure above and below flight level. Flight-level IWC from IKP2 and RASTA agree quite well (top, left), but RASTA often reveals higher IWC below flight level (top and bottom right). Very dense clouds may attenuate the radar beam making IWC retrieval i m p o s s i b l e , a n d e ve n higher IWC may exist b e n e a t h t h e l e v e l o f attenuation.BTD is somewhat indicative of HIWC but remains constant through much of this flight.HIWC was often encountered in the vicinity of detected OTs.High COD is also closely associated with HIWC, but reliable retrievals are only available during the day. Ice particle effective radius exhibits n o u s e f u l s i g n a l f o r H I W C identification. from GOES-13) of (a) visible reflectance, (b) BTW, (c) BTD (d) COD, (e) ΔBTW, (f) visible texture rating (magenta contours) and OT probability rating (color shading), (g), distance to nearest OT, and (h) PHIWC for Cayenne flight #15 on 16 May, 2015." CloudSat overpass of Cayenne flight domain on May 16, 2015 at 1720 UTC. IWC retrievals from the CloudSat 2B-CWC-RO product are shown here (top panel). Flight-level and column-maximum IWC were extracted and plotted along with in-situ TWC and PHIWC (bottom panel). ! Image time and aircraft position on the figure at left correspond to boxed time period (0019 UTC) on time series (right).

Acknowledgements

This research was supported by the NASA High Ice Water Content Program. TWC and RASTA datasets are distributed by the NCAR/UCAR Earth Observing Laboratory. CloudSat Standard Data Products are distributed by the CloudSat Data Processing Center located at the Cooperative Institute for Research in the Atmosphere (CIRA) at Colorado State University in Fort Collins.

PHIWC Validation/Verification Methods

Given that "HIWC" is a loosely-defined term, validation of HIWC diagnostic tools is a challenge. "Impulse events" can be minimized by avoiding strong updrafts as previously shown, but exposure of aircraft engines to low or moderate IWC for a long duration is also an important scenario to consider. PHIWC > 0.6 was obtained for most TWC > 1 g m -3 using the fuzzy logic method presented here, but many moderate TWC were also given similar values. ~90% of TWC < 0.25 had PHIWC < 0.6 (top).

Computation of HIWC probability of detection (POD) and false-alarm rates (FAR) are one method for evaluating performance. Receiver-operator characteristic (ROC) curves are a convenient method to display both quantities. POD and FAR are shown here as functions of different PHIWC thresholds used to signify HIWC. Lower thresholds result in higher POD but also higher FAR.

During the Cayenne campaign there was a coordinated flight with a CloudSat overpass at ~1720 UTC on May 16, 2015. This provided a unique opportunity to compare in-situ measurements with satellite retrievals and to validate HIWC diagnostic products.


J.-A Désidéri

R Duvigneau
email: regis.duvigneau@inria.fr

Parametric optimization of pulsating jets in unsteady flow by Multiple-Gradient Descent Algorithm (MGDA)

Keywords: Active-flow control, time-dependent Navier-Stokes equations, finitevolume schemes, sensitivity equations, multi-objective differentiable optimization, descent methods, robust design 1

Two numerical methodologies are combined to optimize six design characteristics of a system of pulsating jets acting on a laminar boundary layer governed by the compressible Navier-Stokes equations in a time-periodic regime. The flow is simulated by second-order in time and space finite-volumes, and the simulation provides the drag as a function of time. Simultaneously, the sensitivity equations, obtained by differentiating the governing equations w.r.t. the six parameters are also marched in time, and this provides the six-component parametric gradient of drag. When the periodic regime is reached numerically, one thus disposes of an objective-function, drag, to be minimized, and its parametric gradient, at all times of a period. Second, the parametric optimization is conducted as a multi-point problem by the Multiple-Gradient Descent Algorithm (MGDA) which permits to reduce the objective-function at all times simultaneously, and not simply in the sense of a weighted average.

fects at a moderate energy expense. It has been a growing research area for the last decades, since this approach demonstrated its ability to improve aerodynamic performance [START_REF] El-Hack | Flow control: fundamentals and practices[END_REF], for a large range of applications. It is especially appealing in case of separated flows, for which natural instability phenomena can be efficiently exploited to manipulate flow characteristics using periodic flow excitation.

In this context, a major difficulty is related to the choice of actuation parameters, such as excitation frequency, amplitude, location, to obtain the expected flow response. In cases implying a single isolated actuator, it is relatively easy to carry out an experimental or numerical study to determine efficient control parameters. However, in the perspective of industrial applications involving hundreds of actuators, this task is far from being straightforward and the use of an automated optimization strategy is thus proposed, in the spirit of previous works [START_REF] Duvigneau | Optimal location of a synthetic jet on an airfoil for stall control[END_REF][START_REF] Duvigneau | Comparison of turbulence closures for optimized active control[END_REF][START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF].

The application of an optimization procedure to such problems is faced to the following difficulties: first, the choice of the optimization algorithm is conditioned by the huge computational time of the unsteady-flow simulation, and second, it is necessary to consider several objectives concurrently. Typically, the improvement of the single time-averaged performance is usually not satisfactory for realistic applications. Secondly, sensitivity analysis is tedious in the context of unsteady flows, due to the backward integration of the adjoint equation, which requires the storage, or partial storage / partial re-computation, of the unsteady solution.

The proposed work is based on two methodological ingredients to overcome the difficulties described above: the Sensitivity Equation Method (SEM) for unsteady flows on one side, which allows to compute the gradient of a cost-functional with respect to (w.r.t.) control parameters at any time using a forward time-integration, and the Multiple Gradient Descent Algorithm (MGDA) on the other side, which is an extension of the classical steepest-descent method to multiobjective problems and permits to compute a descent direction common to a possibly-large set of costfunctions.

Problem description: optimization of pulsating jets

We consider as model problem the two dimensional compressible flow over a flat plate equipped with three periodically oscillating jets (see Fig. 1). The Reynolds number based on the length h is R = 10 3 and the flow is laminar, while the Mach number is M = 10 -1 . For the three jets, the crosswise velocity is imposed as:

v k (x,t) = A k sin(2πN k t + ϕ k )ζ (x) k = 1, 2, 3, (1)
where ζ (x) corresponds to a squared sine distribution. The jet frequencies are set to the fixed values

N 1 = N ∞ , N 2 = 2N ∞ and N 3 = 1/2N ∞ , with N ∞ = u ∞ /h.
x = {A 1 , A 2 , A 3 , ϕ 1 , ϕ 2 , ϕ 3 }.
Initial values are chosen somewhat arbitrarily as x 0 = {u ∞ , 3/2 u ∞ , 2 u ∞ , 0, π/4, 3π/4}. The grid employed in this study counts 111 161 nodes (see Fig. 2). The initial solution corresponds to uniform flow variables based on inlet conditions. The time step is set to ∆t = 1/(400N ∞ ). The unsteady flow tends rapidly to a periodic regime, whose period corresponds to the lowest actuation frequency 1/2N ∞ and is thus described by 800 time-steps. As transient effects have vanished, one period is defined as observation interval (see Fig. 3). Instantaneous velocity fields are shown in Figs. 4-5 as illustration.

We consider a set of objective-functions { f j (x)} j=1,...,m , defined as the values of the drag J , estimated at discrete times {t j } j=1,...,m , chosen in the observation interval:

f j (x) = J (W j ) with W j = W(x,t j ) ∀ j ∈ {1, • • • , m}, (2)
where x ∈ R n represents the vector of the n = 6 control parameters and W the flow variables. The objective of this work is to reduce simultaneously these m costfunctions. The following sections describe how the gradient ∇ x f j of f j w.r.t. x is evaluated, and the optimization algorithm proposed to conduct the simultaneous optimization of these functions. 3 Sensitivity analysis for unsteady flow

Method

The governing flow equations are written in conservative form as follows:

∂ W ∂t + ∇ • F = ∇ • G , (3)
where W = (ρ, ρu, ρv, ρe) is the vector of conservative mean-flow variables, ρ is density, u and v are the velocity components, and e the total energy per unit mass;

F = (F x (W), F y (W)) and G = (G x (W), G y (W)
) are the vectors of convective and diffusive fluxes respectively. Here ∇ stands for the gradient w.r.t. the spatial Cartesian coordinates x and y, and (∇.) for the divergence operator. The pressure p is obtained from the perfect-gas state equation:

p = ρ(γ -1)(e - u 2 + v 2 2 ) = ρ(γ -1)e i (4)
where γ = 7 5 is the ratio of the specific heats for diatomic gas, and e i the internal energy.

The inviscid fluxes are given by:

F x (W) =     ρu ρu 2 + p ρuv ρu(e + p ρ )     F y (W) =     ρv ρvu ρv 2 + p ρv(e + p ρ )     . (5)
The viscous fluxes are written as:

G x (W) =     0 τ xx τ yx uτ xx + vτ yx -q x     G y (W) =     0 τ xy τ yy uτ xy + vτ yy -q y     , (6)
where τ is the symmetric viscous stress tensor and q the heat flux.

We can now introduce the sensitivity field W , which is defined as the derivative of the flow solution W w.r.t. a given control parameter a, component of x:

W = ∂ W ∂ a . ( 7
)
The equations governing the sensitivity field can be obtained by differentiating (3) w.r.t. a:

∂ ∂ a ∂ W ∂t + ∂ ∂ a (∇ • F ) = ∂ ∂ a (∇ • G ) . (8)
By switching the derivatives w.r.t. a and those w.r.t time or space coordinates, one obtains:

∂ ∂t ∂ W ∂ a + ∇ • ∂ F ∂ a = ∇ • ∂ G ∂ a , (9)
or:

∂ W ∂t + ∇ • F = ∇ • G , (10)
which is formally similar to (3), by introducing the sensitivity of the convective flux F = (F x (W, W ), F y (W, W )) and the sensitivity of the diffusive flux

G = (G x (W, W ), G y (W, W )).
The sensitivity of the convective fluxes can be expressed as:

F x (W, W ) =     (ρu) (ρu) u + (ρu)u + p (ρu) v + (ρu)v (ρu) (e + p ρ ) + (ρu)(e + ( p ρ ) )     (11)
F y (W, W ) =     (ρv) (ρv) u + (ρv)u (ρv) v + (ρv)v + p (ρv) (e + p ρ ) + (ρv)(e + ( p ρ ) )     . ( 12
)
The sensitivity of the diffusive fluxes reads:

G x (W, W ) =     0 τ xx τ yx u τ xx + v τ yx + uτ xx + vτ yx -q x     (13) G y (W, W ) =     0 τ xy τ yy u τ xy + v τ yy + uτ xy + vτ yy -q y     , ( 14
)
where τ is the sensitivity of the viscous stress tensor and q the sensitivity of the heat flux. The boundary conditions for the sensitivity equations are obtained by differentiating the boundary conditions applied to the flow.

Since the flow and sensitivity equations are formally similar, both are solved using the same finite-volume approach [START_REF] Duvigneau | A sensitivity equation method for unsteady compressible flows: Implementation and verification[END_REF], based on a second-order vertex-centered discretization scheme. Temporal integration relies on a second-order implicit backward method, with a dual time-stepping technique. Note that the implicit part of the scheme is the same for the flow and sensitivity equations, since both involve the same Jacobian matrix.

The sensitivity equation depends on the parameter a, component of x of interest. Therefore, six sensitivity equations have to be solved, possibly in parallel, to estimate the components of the gradient for all m cost-functions:

∇ a f j = ∂ W J (W j ) • W (t j ) ∀ j ∈ {1, • • • , m}. (15)
We emphasize that if m is large, the solution of six sensitivity equations is far more cost-efficient than solving m adjoint equations backward in time. Additionally, the memory-storage requirement remains moderate.

Verification

To verify the implementation of the sensitivity equations, we compute a neighboring solution according to a first-order extrapolation W(a) + W δ a and compare it with the solution W(a + δ a). This exercise is conducted in a simplified case including a single jet, a being the jet amplitude A 1 . Figs. 6-7 provide illustrations for the flow fields at selected times and Fig. 8 for the resulting drag history, for a perturbation of the jet amplitude δ A 1 = A 1 /4. A similar exercise has been achieved for the phase, to fully verify the gradient estimation. Fig. 8 Linear extrapolation of the drag w.r.t. jet amplitude A 1 : reference drag J (A 1 ) in green, extrapolated drag

J (A 1 ) + J δ A 1 in red, non-linear perturbed drag J (A 1 + δ A 1 ) in blue, for δ A 1 = A 1 /4.

Multiobjective descent algorithm MGDA

Equipped with procedures for calculating the objective-functions and their gradients, we now turn to the issue of constructing the multi-objective optimization method.

The Multiple-Gradient Descent Algorithm (MGDA) was originally introduced in [1] and [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF] to solve general multi-objective optimization problems involving differentiable cost-functions. Variants were proposed in [START_REF] Désidéri | Numerical Methods for Differential Equations, Optimization, and Technological Problems[END_REF], but more recently the algorithm was slightly revised in [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF] to apply to cases where the number m of objectivefunctions exceeds the dimension n of the working design space. We recall here the basic definition of the revised version and provide some details about the application to the present parametric optimization.

Multi-Objective problem statement

Let m and n be two arbitrary integers, and consider the multi-objective optimization problem consisting in minimizing m differentiable objective-functions { f j (x)} in some open admissible domain Ω a ⊆ R n ( j = 1, . . . , m; f j ∈ C 1 (Ω a )). Given a starting point x 0 ∈ Ω a and a vector d ∈ R n , one forms the directional derivatives

f j = [∇ x f j (x 0 )] t d (16)
where ∇ x is the symbol for the gradient w.r.t. x and the superscript t stands for transposition. One seeks for a vector d such that

f j > 0 (∀ j) . ( 17
)
If such a vector d exists, the direction of vector (-d) is said to be a local descent direction common to all objective-functions. Then evidently, infinitely-many other such directions also exist, and our algorithm permits to identify at least one.

Convex hull, two lemmas and basic MGDA

We recall the following :

Definition 1. The convex hull of a family of m vectors {u j } ( j = 1, . . . , m; u j ∈ R n ),
is the set of all their convex combinations:

U = u ∈ R n such that u = m ∑ j=1 α j u j ; α j ∈ R+; m ∑ j=1 α j = 1 . (18)
Then, we have :

Lemma 1. Given an n×n real-symmetric positive-definite matrix A n , the associated scalar product

u, v = u t A n v (u, v ∈ R n ), (19)
and Euclidean norm

u = u t A n u, (20)
the convex hull U admits a unique element ω of minimum norm.

Proof. -Existence : U is closed and . is a continuous function.

-Uniqueness : suppose that ω 1 and ω 2 are two realizations of the minimum µ = arg min u∈U u so that µ = ω 1 = ω 2 and let

ω s = 1 2 (ω 2 + ω 1 ) , ω d = 1 2 (ω 2 -ω 1 )
, so that:

ω s , ω d ) = 1 4 ω 2 + ω 1 , ω 2 -ω 1 = 1 4 ω 2 2 -ω 1 2 = 0 .
Hence ω s ⊥ ω d , and since ω s ∈ U, ω s ≥ µ, and :

µ 2 = ω 2 2 = ω s + ω d 2 = ω s 2 + ω d 2 ≥ µ 2 + ω d 2 =⇒ ω d = 0 .
Lemma 2. The minimum-norm element ω defined in Lemma 1 satisfies :

∀u ∈ U, u, ω ≥ ω 2 . (21)
Proof. Let u ∈ U, arbitrary. Let δ = uω; by convexity of U :

∀ε ∈ [0, 1], (1 -ε)ω + εu = ω + εδ ∈ U ,
and by definition of ω, ω + εδ ≥ ω , that is :

ω + εδ , ω + εδ -ω, ω = 2ε ω, δ ) + ε 2 δ 2 ≥ 0 ,
and this requires that the coefficient of ε be non-negative.

Then consider the case where

u j = ∇ x f j (x 0 ) (∀ j). (22)
If the vector ω defined in Lemma 1 is nonzero, the vector

d = A n ω (23)
is also nonzero, and is a solution to the problem stated in ( 16)-(17) since by virtue of Lemma 2:

u j , ω = u t j A n ω = u t j d ≥ ω 2 > 0. ( 24
)
The situation in which ω = 0, or equivalently,

∃α = {α j } ∈ R + m such that m ∑ j=1 α j ∇ f j (x 0 ) = 0 and m ∑ j=1 α j = 1, (25)
is said to be one of "Pareto-stationarity". The relationship between Pareto-optimality and Pareto-stationarity was made precise by the following [3]- [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF] Theorem 1. If the objective-functions are differentiable and convex in some open ball B ⊆ Ω a about x 0 , and if x 0 is Pareto-optimal, then the Pareto-stationarity condition is satisfied at x 0 .

Hence, the Pareto-stationarity condition generalizes to the multi-objective context, the classical stationarity condition expressing that an unconstrained differentiable function is extremal.

We now return to the non-trivial case of a point x 0 that is not Pareto-stationary and we suppose that the vectors ω and d (ω = 0; d = 0) have been identified (see next subsection). Then we define MGDA as the iteration which transforms x 0 in

x 1 = x 0 -ρd (26)
where ρ > 0 is some appropriate step-size. Thus MGDA is an extension to the multiobjective context of the classical steepest-descent method, in which the direction of search is taken to be the vector d defined above. At convergence, the limiting point is Pareto-stationary.

We now examine how can the vector d be computed in practice.

QP formulation and hierarchical Gram-Schmidt orthogonalization

By letting

ω = m ∑ j=1 α j u j = Uα ( 27
)
where u j = ∇ f j (x 0 ), U is the n × m matrix whose jth column contains the n components of vector u j , the identification of vector ω can be made by solving the following Quadratic-Programming (QP) problem for the unknown vector of coefficients α = {α j }:

ω = arg min α∈R m 1 2 α t Hα (28)
subject to:

α j ≥ 0 (∀ j), m ∑ j=1 α j = 1, (29)
where H = U t A n U. Note that if vector ω is unique, vector α may not be.

If the family of gradients is linearly-independent, which requires in particular that m ≤ n, it is possible to choose the scalar product, through the definition of matrix A n , in such a way that these gradients form an orthogonal basis of their span. Then vector ω is explicitly determined by the orthogonal projection of 0 onto the convex hull U:

α j = 1 u j 2 ∑ m k=1 1 u k 2 (30)
In the inverse case where m > n (and even m n), focus of interest presently, let r ≤ n < m be the rank of the family of gradients. Using first the standard Euclidean scalar product (A n = I n ), the Gram-Schmidt orthogonalization process stops in r steps and produces an orthogonal basis (in the usual sense) {v j } ( j = 1, . . . , r), that span a subspace G of dimension r. The orthogonal vectors, {v j }, are calculated from a subfamily of the original vectors, {u j }, j ∈ J, where J is a subfamily of r indices from 1 to m. In this process, a hierarchical principle was introduced in [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF] to select these indices in such a way that the cone bounded by the reduced family {u j } ( j ∈ J) be as large as possible to contain the directions, in the most favorable situation, of all the other gradients, unused in the Gram-Schmidt process by redundancy. When this occurs, the subfamily {u j } ( j ∈ J) not only is a basis of the subspace G , but also, its convex hull contains all the directions of interest. In the more general case where the directions of some gradients among the unused vectors {u j } ( j / ∈ J) are not in the cone, we resort to the QP-formulation, but expressed after changing the basis to become {u j } ( j ∈ J) and by choosing a new scalar product to make this subfamily orthogonal. The technical steps are the following [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF]:

• Once the n × m matrix U is formed with the components of the given gradients {u j } (u j ∈ R n , j = 1, . . . , m), these gradients, or column-vectors are made (physically) dimensionless, by component-wise normalization to form the initial matrix G : ∀i ∈ 1, . . . , n :

s i = max j u i, j , S = Diag s i , G = S -1 U. ( 31
)
This normalization is essential to make the subsequent calculations of scalar products physically meaningful and computationally well-balanced. • Throughout the Gram-Schmidt process, columns of the G matrix are permuted by the hierarchical selection of basis vectors. Upon exit, the actually used gradients are placed in the first r columns. The corresponding n × r leftmost block of the final matrix G is then denoted G. Note that in the present version of the Gram-Schmidt process, the computed orthogonal vectors {v j } are not normalized to unity, but in a special way for which r directional derivatives are equal [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF]. Let these vectors be stored in matrix V, and define the following diagonal matrix:

∆ = Diag v t j v j . (32)
Then:

A n = W t W + (I -Π ) 2 , W = G t G -1 G t V∆ -1 V t , (33)
where Π = V∆ -1 V t is the projection matrix onto subspace G [START_REF] Désidéri | Révision de l'algorithme de descente à gradients multiples (MGDA) par orthogonalisation hiérarchique[END_REF].

In this way, the QP-formulation is well-conditioned and easily solved by a library procedure. We have used the procedure qpsolve from the Scilab library which is equivalent to the quadprog procedure from the MATLAB library. As a result of this, exactly r directional derivatives { f j } are equal, and by experience, the remaining ones differ only slightly.

Results

For each control parameter, component of x = {A 1 , A 2 , A 3 , ϕ 1 , ϕ 2 , ϕ 3 }, the SEM is applied to obtain sensitivity fields. We provide as illustration (see Figs. We aim now at determining the vector of parameters x = {A 1 , A 2 , A 3 , ϕ 1 , ϕ 2 , ϕ 3 } that reduces simultaneously the 800 cost-functions associated with the observation period. To reduce somewhat the computational complexity without altering greatly the transient behavior, the MGDA approach is applied to only m = 20 homogenized gradients, obtained by averaging the gradients by time-intervals of 40 time-steps. As a result, one disposes of a common descent direction associated with vector d satisfying (17).

Once the vector d is determined, a practical step-size ρ must be estimated. For this, we first note that a natural scale for the variations of a time-dependent objectivefunction is given by its standard deviation, σ , a more significant value than its average which can be 0. Then if δ x = -ρd, the variation of the objective-function average can be estimated as -ρ ḡ.d where ḡ is the average gradient. Thus, a meaningful reference step-size can be defined by the condition:

ρ ḡ.d = σ (34) where σ = 1 m ∑ m j=1 f j -f 2 , f = 1 m ∑ m j=1 f j , and ḡ = 1 m ∑ m j=1 ∇ f j (x 0 ). This gives: ρ = σ ḡ.d . (35)
In practice, in the present experiments, we have used the step-size ρ = 1 10 ρ to update the control parameters by the descent method, and this resulted in a successful iteration, stable and effective.

The history of the drag in the observation period is represented in Fig. 12, from the baseline flow to full convergence of the optimization approach. As expected, each update of the design vector has resulted in a diminished drag over the entire observation period. As it can be noticed, some points in time are more critical than others. In contrast, when one applies, more classically, the steepest-descent method to the time-averaged cost-function J , by setting the search direction to the average gradient, an increase of the drag can be observed at some times, as illustrated in Fig. 13. Finally, also note that actuation permits a significant drag reduction w.r.t. the case without suction/blowing for which the value of drag is indicated on the figure by a dotted horizontal line.

Finally, a second exercise is conducted: the drag values computed over the last 40% of the observation period only are considered as optimization criteria in the MGDA approach (in this interval, the drag is especially high for the baseline flow). The history of the drag in the observation period is represented in Fig. 14, for the first 16 iterations of the optimization algorithm. As expected, a more significant decrease is achieved during the last 40% of the observation period, whereas the drag is free to vary in the first 60%, and in fact, increases.

Conclusion

In this work, we solved for demonstration an exercise of active-flow control in which drag over a flat plate has been reduced by three pulsating jets acting on the boundary layer. The flow, governed by the time-dependent compressible Navier-Stokes equations, has been simulated numerically by second-order in time and space finitevolumes, yielding, when the periodic regime is achieved, drag as a function of time.

The simultaneous solution of the sensitivity equations has provided additionally the six-component gradient of drag w.r.t. the design characteristics of the jets. The accuracy of these gradients has been verified by comparison with finite-differences via fine-mesh computations. By the Multiple-Gradient Descent Algorithm (MGDA), a direction of search has been identified permitting to reduce drag at all times of the period, or a selected segment of it. The process was repeated iteratively, and at all intermediate steps of the optimization process as well as at convergence, the drag was effectively reduced over the entire observation time-interval.

Hence, we dispose of a numerical optimization tool whose efficacy is demonstrated uniformly, that is, over a possibly-large range of operational conditions, here different discretization times. This contrasts with more classical approaches in which a single functional, usually defined as a somewhat arbitrary weighted average, is minimized at the risk of a degradation of certain elements composing the average. This method is currently being extended to solve more general robust design problems.

Fig. 1

 1 Fig. 1 Problem description.

Fig. 2

 2 Fig. 2 Computational mesh.

Fig.

  Fig. History of drag for initial parameters and observation area.

Fig.

  Fig. Snapshot of the streamwise velocity field.

Fig.

  Fig. Snapshot of the crosswise velocity field.

Fig. 6

 6 Fig. 6 Linear extrapolation of streamwise velocity field u w.r.t. jet amplitude A 1 for blowing (top) and suction (bottom) phases: reference state u(A 1 ) in green, extrapolated state u(A 1 ) + u δ A 1 in red, non-linear perturbed state u(A 1 + δ A 1 ) in blue, for δ A 1 = A 1 /4.

Fig. 7

 7 Fig. 7 Linear extrapolation of pressure field p w.r.t. jet amplitude A 1 for blowing (top) and suction (bottom) phases: reference state p(A 1 ) in green, extrapolated state p(A 1 )+ p δ A 1 in red, non-linear perturbed state p(A 1 + δ A 1 ) in blue, for δ A 1 = A 1 /4.

  9-10) instantaneous sensitivity fields of the velocity w.r.t. the amplitude for the first jet. The derivatives of the drag w.r.t. the control parameters are then computed at all timesteps, yielding 800 values of cost-function and gradient for the whole observation period, as illustrated in Fig.11.

Fig. 9

 9 Fig. 9 Sensitivity of streamwise velocity w.r.t. first jet amplitude.

Fig. 10 Fig. 11

 1011 Fig. 10 Sensitivity of crosswise velocity w.r.t. first jet amplitude.

  The jet amplitudes and phases are considered as control parameters

  Evolution of the drag history w.r.t. optimization iterations: case of a mean direction descent (blue: initial, red: after 2 iterations, black: intermediate iteration).Fig.14Evolution of the drag history w.r.t. optimization iterations: case of MGDA approach based on the last 40% of the observation period (blue: initial, red: after 7 iterations, black: intermediate iterations).

	0.07
		initial
		Mean direction 2 it
	0.06
	0.05
	drag
	0.04
	0.03
	0.02
	13.5	14	14.5	15
		time
	0.03 0.04 0.05 0.06 0.07 0.03 0.04 0.05 0.06 Fig. 13 13.5 drag drag	initial MGDA (full interval) converged without blowing-suction 14 initial mgda (40% interval) it 16 without blowing / suction	14.5	15
	0.02	time
	13.5	14	14.5	15
		time
	Fig. 12 Evolution of the drag history w.r.t. optimization iterations: case of MGDA approach (blue:
	initial, red: final, black: intermediate iterations).