domain: spi phys

text:
Ghislain

Pierre B Lokman

Daniel Tobias

Fabien Julie Olivier

Adrien Alexandre

Bastien

Benjamin

Charlotte

Alexandros

Maxime, Thibault, Hala, Virgile, Julien, Tanguy, . . . Rémi

Pour commencer, je remercie chaleureusement les membres du jury qui ont participé à ma soutenance. J'ai apprécié vos questions, remarques et commentaires bienveillants qui ont alimenté la discussion tout au long de la soutenance. Une mention spéciale pour les trois rapporteurs Eric Laurendeau, Pierre Lubin et Christian Ruyer-Quil qui ont accepté la lourde tâche d'être rapporteur pour mon manuscrit. Je vous suis reconnaissant d'avoir consacré du temps pour examiner mes travaux.

Je remercie également le directeur du DMPE, Pierre Millan, qui a facilité mon retour en 2012 sur Toulouse pour travailler sur les problématiques liées au givrage. Jamais je n'aurais imaginé à l'époque que ce serait un virage décisif pour moi. Un grand merci également à Pierre Gajan, mon chef d'unité, qui a toujours favorisé l'épanouissement scientifique chez les jeunes de son groupe.

Evidemment, un grand merci à Philippe Villedieu, le chef d'orchestre de l'équipe givrage. Ta curiosité scientifique est communicative et ta rigueur mathématique rassurante ! Je remercie également tous les acteurs du givrage à l'ONERA : Emmanuel,

Chapter 1

Introduction

This manuscript, written in order to obtain the "Habilitation à Diriger des Recherches" (HDR), presents the works I have been carrying out at ONERA since 2009 as a research engineer. It reflects the two steps in my career at ONERA. First of all at the DSNA department (now DAAA) in Châtillon where my topic was optimization and gradient computation by an adjoint method. It was also the opportunity for me to experience high order numerical methods like isogeometric analysis. The real watershed for me was my transfer to the DMAE department (now DMPE) headed by P. Millan in Toulouse. At that time, icing activities had just been transferred to the DMAE under the supervision of P. Villedieu. The new team I joined from the beginning and which has grown over the years, had to address requirements inforced by the emergence of new certification rules such as SLD or ice crystal icing. Supported by French industry and European funding, icing activities have enjoyed significant dynamism at the national and international level over the past decade. Thus, the requirements for new and ever more reliable databases as well as modern computational icing suites have been strong. It is in this context that ONERA had to investigate many research subjects related to icing, such as complex wall impacts (SLD, ice crystals), wall roughness characterization, thermodynamic balance, thin film runback, ice protection systems, hydro/ice-phobic materials, mechanical characterization of accreted ice, advanced numerical methods like remeshing techniques, . . . ONERA has even invested in the building of a research icing wind tunnel. I have been involved in most of these topics, either directly through the coding of the 2D icing suite IGLOO2D and the derivation of the models integrated in the code, or through the supervision of interns and PhD students.

The manuscript is organized as follows. It is composed of four parts. The two first ones are dedicated to my icing activities at the DMPE. More precisely, the first part is focused on the socalled cold icing conditions with supercooled droplets and ice crystals. The second part deals with the modeling of electro-thermal ice protection systems. The third one briefly summarizes my PhD results and then is devoted to my works at the DSNA (adjoint method and isogeometric analysis). Each part is itself divided into chapters dedicated to specific topics. Each of these chapters is introduced by a few lines in which I list the publications (papers, conference papers and oral presentations) I have written and which are directly related to the topic in question. If the reported results are obtained within the framework of PhD students I supervised, this is emphasised during these few introduction lines. Finally, each chapter is ended by a conclusion in which the perspectives are presented. The general conclusions and especially the general perspectives are detailed in the final chapter of this document. The fourth and last part of the manuscript (Appendices) is dedicated to non-technical information namely a detailed CV, a publication list and five relevant papers. 9

Part I

Ice accretion: basic concepts and modeling

This part is composed of four chapters:

• Chapter 2: a brief introduction to icing to clarify the scientific context of this activity. It should be mentioned that it also introduces concepts presented in the second part (Sec. II).

• Chapter 3: icing conditions with supercooled liquid water droplets. Emphasis is put on large droplets (SLD) and their interactions with a wall.

• Chapter 4: ice crystal icing conditions.

• Chapter 5: description of the 2D icing suite IGLOO2D

Chapter 2

A brief introduction to icing

Why a specific introduction for my icing-related activities ? To explain how icing is at the crossroads among several scientific domains. Icing is actually not a fundamental subject that could be studied at university, such as thermodynamics or fluid mechanics, for instance. Moreover, icing is often a source of concern raised by aircraft or engine manufacturers themselves (like air safety concern or ice protection system improvement). Therefore, the related challenges may appear to be, at first glance, only of industrial nature with technical and technological solutions. However, the accurate computation of an ice shape or the estimated energy required from an ice protection system to prevent its growth are actually scientific challenges. Among the key topics related to icing, one can highlight:

Thermodynamics

The starting point for ice build-up involves ice particles and supercooled liquid water droplets which exist in liquid form at temperatures below 0 • C. Supercooling requires water to be pure and free of nucleation sites. An abrupt crystallization of these droplets may happen when they hit a flying aircraft. This change from liquid to solid state is associated with a release of latent heat. Ice accretion is driven by two conservative equations. Firstly, as far as mass is concerned, a balance is reached between the part of the water that freezes, the one that remains in liquid state and the one that evaporates. Regarding energy conservation, a competition occurs between the source terms (release of latent heat by the impact of supercooled water droplets, heating provided by an ice protection system, . . . ) and the sink terms (energy extracted from the substrate by the sticking of fully frozen ice particles, convective heat transfer between the iced wall and the cold air boundary layer, evaporation and sublimation,. . . ). The balances as well as the magnitudes for each of the contributing terms are detailed in Sec. 3.1.

Dispersed phase

Before contributing to the thermodynamic balance on the aircraft, engine or probe walls, the droplets and ice particles are conveyed into the gaseous flow. Two particular phenomena can then be focused on.

Firstly, the particle transportation step itself. Regarding droplets with diameters larger than 50µm (see Sec. 3.2 on SLDs), the large air velocity gradients observed near the stagnation point of an airfoil may cause the droplets to atomize due to the local increase in the value of the Weber number. Regarding ice crystals (Sec. 4), the aerodynamic drag and lift coefficients of the particle must be 15 adapted to take into account its irregular shape which cannot be assumed to be spherical. In addition, due to the melting phenomenon that occurs when an iced particle enters a hot environment, its shape may change over time.

Complex particle-wall interactions are also an active research topic. Although this has been widely investigated for droplets at low impact velocity, many physical phenomena still remain poorly understood for the typical configurations encountered in aeronautics (impact velocities ∼ 100 m/s). Hot points are mass losses (splashing) due to partial deposition on the wall, secondary (re-emitted) droplet characterisation (diameter, direction and velocity) and thresholds among splashing, spreading and bouncing regimes. Obviously, the surface property of the wall (roughness, substrate composition, wall temperature, liquid water available on the wall, . . . ) is of paramount concern in determining the nature of the impact. The same challenges can be identified with ice particles and the additional solid phase to be considered.

Aerodynamics and heat transfer

Computing the aerodynamic flow field around a frozen airfoil can be tricky, especially in the case of ice shapes with large horns and boundary layer separation. One of the objectives is to estimate the convective heat transfer coefficient for the computation of the wall heat flux as a main sink term for the thermodynamic balance (see above). In the presence of accreted ice, the boundary layer grows on a rough wall. The challenges are of two kinds. First of all, it is necessary to model the influence of roughness on the computation of the wall friction coefficients C f and the heat transfer coefficient h t . These models, generally well designed for smooth surfaces, must be adapted for rough walls. Therefore, and this is the second difficulty to be addressed, local roughness due to the accreted ice must be characterized from a reduced number of global parameters as inputs for the models for C f and h t . The similitude between real roughness (i.e. the one observed with accreted ice) and academic roughness (i.e. the one used as input data for the models) is based on the equivalent sand grain thickness k s . The definition of k s and more generally the roughness influence on the boundary layer still remains a very active research field.

Another point is the development of a boundary layer above a heated liquid film. The resulting high evaporation rates can locally saturate the boundary layer with steam, which retroactively influences the evaporation rates themselves. This coupling between the boundary layer and the evaporating liquid film is an additional modeling difficulty.

Mechanics and materials

Accreted ice is a fragile material whose physical properties are poorly known. These depend on many parameters, including ice temperature, liquid water content and porosity. Experimental studies for the characterization of laboratory formed ice are available. They are more scarce for ice obtained directly by accretion on aerodynamic bodies.

Regarding ice shedding (Sec. 6.4), in addition to the intrinsic ice properties necessary to model cohesive crack growth (bulk failure), the ice/substrate interaction properties are required to characterize crack growth in adhesive interfacial debonding. Wall properties are therefore additional parameters to be taken into account.

Regarding ice crystal icing, the poor estimation of the ice properties is one of the major hurdles for the modeling of the erosion phenomenon (Sec. 4.4.2). For instance, liquid water content of the slushy ice deposits is not yet predictable, which is annoying since it is one of the relevant parameters to characterize erosion.

Wetting/coating

In the framework of ice protection systems, hydrophobic or icephobic coatings may be preferred to prevent ice build-up. Wetting and coating are both related to the modeling of the interaction between a liquid water film or an ice layer with the wall. A typical application is the transition of a continuous liquid film into rivulets or isolated droplets (Sec. 7). Indeed, predicting the wet to dry surface ratio is a key input for estimating the heat power required to make a liquid film evaporate or prevent it from freezing further. To do this, the contact angle θ s is the main input parameter to describe the contact line at the interface among the gas, the liquid film and the wall. However, θ s itself depends on the local properties of the wall such as cleanliness and roughness. The heterogeneity of θ s on the wall allows different forward and backward angles for an isolated droplet (hysteresis phenomenon), which explains the shape of a sessile droplet while facing air shear flow.

Passive ice protection systems, with the use of coatings, are a promising way forward (Sec. 6).

It is within this scientific framework that my main research activity has been focused on icing modeling since my transfer to the DMAE (now DMPE) department at the end of 2011 (Fig. 2.1). Chronologically, this activity corresponds to the second part of my career. My objectives have been:

• To develop models able to address the requirements inforced by the emergence of new certification rules (Appendix O for Supercooled Large Droplets and Appendix D for mixed phase and ice crystal icing).

• To develop a new 2D icing tool (called IGLOO2D) based on a modular advanced software architecture.

• To address new challenges such as ice protection system modeling with runback liquid films and rivulet formation.

The scientific approach of this work is in line with the mission of ONERA, namely the development of models and methods for aerospace industrial partners. This pragmatic need has sometimes led me to focus on models of a correlating nature. However, enhancing the underlying scientific basis of these models has always been the main concern of my studies.

Chapter 3

Supercooled water icing As introduced in Sec. 2, supercooled water icing is observed with supercooled liquid water droplets which exist in liquid form at temperatures below 0 • C and which freeze when they hit a wall. This chapter is divided into two main sections. Firstly, in Sec. 3.1, focus is on small droplets with a diameter of less than 50 µm. This is the opportunity to introduce the mass and energy balances which are basic references for the other chapters to describe ice buildup. Secondly, in Sec. 3.2, emphasis is on supercooled large droplets (SLD) with a diameter larger than 50 µm. The main difference with smaller droplets lies in the complex wall/particle interaction with splashing and re-emitted droplets. Derivation of models for SLDs has been the subject of several papers [1,2,3,4,5,6] (see Sec. 3.2 for more details).

From supercooled water to ice: basic formulation and application to small droplets

Ice accretion with small supercooled water droplets (diameters less than 50 µm) is presented here. This corresponds to Appendix C for certification rules [7] where the droplets are supposed to stick entirely to the wall. Ice accretion consists in mass and energy balances (Fig. 3.1) among the incoming/outgoing (impacting, running-back or evaporating) mass rates and enthalpies. The main outputs are the resulting regime (rime, glaze or running wet1 ), the estimation of the ice growth rate ṁacc , the accretion surface total liquid mass fraction f l and the wall temperature T w . The different terms of the balances are where m L , U and E k are respectively the liquid water mass, the internal energy and the kinetic energy of the considered system. Note that the previous system (3.1) has to be adapted in the case of ice crystal icing (Sec. 4.4.1). The accretion surface total liquid mass fraction f l is defined as:

f l = ṁrb,out ṁL dep + ṁrb,in -ṁevs (3.2)
Once again the definition of f l has to be adapted in the case of ice crystal icing (Eq. ( 11) from [10]). Equation (3.1) as a function of the wall temperature T w is non-linear, especially via the evaporation/sublimation term ṁevs [9]. A three-step approach, where each of the regimes (rime ice, running wet and glaze ice) are tested, is used to solve Eq. (3.1) and find the three unknowns, namely T w , ṁacc and ṁrb,out [9]. Figure 3.2 shows the Messinger energy balance for a rime ice configuration. The balance is domi- nated by one source term, namely the liquid water deposed mass enthalpy H L dep and two sink terms, namely the convective heat flux Qc (predominant) and the evaporated/sublimated surface mass rate H evs . Thus, the supercooled water droplets, by releasing their latent heat when they solidify, represent the main source term for heat supply in the Messinger's balance. All this heat must be evacuated by evaporation H evs , and mainly by convective transfer Qc (see Fig. 3.2). Qc is computed as follows:

Qc = h t • (T r -T w ) (3.3)
where T r is the recovery temperature and h t the heat transfer coefficient [11]. Hence the need to evaluate precisely h t , which is far from being accurate, especially for rough walls. So, the predominant term Qc in the Messinger balance is one of the most difficult terms to be accurately modeled.

Supercooled Large Droplets (SLD)

Introduction

In 1994, an ATR-72 crashed at Roselawn, Indiana, USA. It has been speculated that this accident was due to SLDs which can impinge beyond the limits of the protection system due to their high inertia. This crash led to a modification of the regulation rules with the definition of the Appendix O [7] which includes freezing drizzle (M V D < 500µm) and freezing rain (M V D > 500µm) icing conditions. To account for these new rules, both experimental and numerical means of compliance need to be modified and validated. To address this problem, several research projects were launched in North America and in Europe (EXTICE European project) during the last decade. Thanks to these works, icing wind tunnels were updated to reproduce SLD conditions and a first set of models have been developed. As far as numerical means of compliance are concerned, the modeling of SLD is very challenging since it requires taking into account several new physical phenomena which do not need to be considered for small droplets (Appendix C conditions). The most important one is related to the so-called "mass loss effect". It corresponds to the fact that large droplet impingement is often accompanied by splashing (Fig. 3.3) with secondary droplet reemission and thus a reduction of the amount of water that would have been deposited otherwise. [12].

Droplet splashing has received great attention in the research community during the last 30 years. From a theoretical/numerical point of view, the spreading of a drop impacting onto a solid surface under isothermal [13,14] and non-isothermal [15,16,17] conditions has been extensively studied. Although understanding of the nucleation phenomenon is a key point in explaining freezing process, few studies have been conducted in the framework of droplet impact [18]. In the case of a supercooled liquid, fast propagation of single dendrites or a cloud of dendrites through the liquid has been observed [19]. See [20] which, in addition to addressing the impact of single supercooled water droplets onto a smooth ice surface, also provides a complete review of academic studies on the impact of SLDs. The above-mentioned academic studies are fundamental for the detailed understanding of basic physics and the definition of characteristic time and space scales dedicated to the impact of SLDs. However, the considered droplet impacts are most often normal with impact velocities in the order of a few m.s -1 . This is far from the specific icing applications with high impact velocities (around 100 m.s -1 ) and a large range for the impact angle (from normal at the leading edge to grazing near the impingement limit). For this reason, dedicated experiments were performed by Papadakis et al. in the NASA Icing Research Tunnel to get a database in representative conditions [21,22,23,24]. These experimental data combined with some ideas issued from existing impingement models (Bai & Gosman [25] for impacts at low velocity, Trujillo et al. [26] Mundo et al. [27,28], Yarin & Weiss [29]) are used by several authors [30,31,32,33,34,35,2,36,37,2,38,39] in order to derive correlative SLD mass loss models for icing applications. Even if these models yield significantly improved collection efficiency predictions when applied to the experiments from Papadakis et al., they suffer from several important weaknesses:

• They are purely empirical models without any theoretical background, except the fact that they rely on some relevant dimensionless numbers introduced in former studies on droplet impingement [27] [29].

• In all the papers [32,39,2,31] except in [38], there is no explanation on how the proposed mass loss model could be explicitely derived from the experimental data from Papadakis et al. It is only shown that the model allows improving the collection efficiency prediction.

• There exist large discrepancies among the existing model predictions despite the fact that they are all based on the same experimental database.

In the context of PHYSICE and PHYSICE2, our objective was to propose a rationally based methodology to derive a more reliable mass loss model [1] than the existing ones. To do this, we focused on:

• The use of experimental databases (Sec. 3.2.2). Once again the experimental database from Papadakis et al. were used to calibrate and validate our model (Sec. 3.2.2). In particular, the model allows to separate clearly the influence of the normal kinetic energy, which is the dominating effect close to the leading edge, from the influence of the angle of incidence which is the most influent parameter close to the impingement limits. Moreover the model can be used for both SLDs and small droplets. In order to extend the droplet velocity range (limited to 78 m.s -1 in the Papadakis et al. database), experiments were proposed by ONERA in the framework of PHYSICE and PHYSICE2 [40,41,42,43] (Sec. 3.2.2). The impingement limit between bouncing and splashing regimes are studied as well as the role of the surface texture of the impacted wall. Moreover, correlations regarding the characteristics (diameter, velocity and direction) of the secondary (re-emitted) particles are derived from the ONERA database (Sec.

3.2.2).

• The use of Direct Numerical Simulations (DNS) of the impact of a single droplet (Sec. 3.2.3). The approach is challenged by the multi-scale nature of the droplet impact as well as the numerical difficulties associated with the front capturing method of an interface and the associated jump conditions. Thus, the systematic use of DNS as a means of building numerical database is still limited. This topic, strongly related to my PhD works (Sec. 8), has led me to supervise the PhD works of T. Xavier.

Revisited Model for SLD impact onto a dry solid surface

Model calibration and validation from the database of Papadakis et al.

This section is mainly derived from [1] where our objective was to proposed a model for the sticking efficiency defined by:

= β dep β imp = ṁdep ṁimp (3.4)
where β dep and β imp are respectively the collection and impingement efficiencies and which are respectively the dimensionless expressions of ṁdep and ṁimp (see Tab. [21,22,23,24] measured β dep thanks to the blotter paper technique. It was based on a water spray containing a known concentration of blue dye. The models were covered with a heavy weight blotter paper and the amount of dye is then measured via reflectance spectroscopy using a CCD camera. Collection efficiencies β dep were determined on different clean and iced airfoils, for several drop size distributions from MVDs of 11µm to 236µm and for several angles of attack from 0 • to 8 • . The airspeed for all cases was 78 m.s -1 . Since β imp cannot be measured, it is not possible to extract the values of the mass loss coefficient from the experimental data. However β imp can be numerically computed, with a quite good level of accuracy, by using the droplet trajectory solver TRAJL2D from the ONERA icing suite IGLOO2D (Sec. 5.2). Combining results from both experimental databases and numerical simulations is the main originality of our paper [1]. It is possible to simultaneously know the characteristics of the impinging droplets (velocity and diameter), the impinging mass flow rate, the deposited mass flow rate and the mass loss coefficient at each point of the test model. Figure 3.4 shows a comparison between experimental results and our SLD model for the sticking efficiency as a function of the angle of incidence θ (or in an equivalent way as a function of the impact angle α, see Fig. 3.4a). Two kinds of profiles are represented: the MS(1)-317 (Fig. 3.4b) and the NACA23012 (Fig. 3.4c) profiles. Even if the results are scattered due to the numerous experimental and numerical uncertainties (droplet cloud homogeneity, droplet size distribution, measurement inaccuracy near the impingement limit, etc.), two impingement regimes may be clearly distinguished for both the suction and the pressure sides as far as SLDs are concerned (M V D > 50 µm):

1. For small angles of incidence (corresponding to large α), the sticking efficiency is an increasing function of |θ| = 90 • -|α|, which tends to zero when |θ| is close to 0 • (impingement limit).

2. For intermediate and large angles of incidence, the sticking efficiency does not strongly depend on θ, and even seems to be almost constant in several cases. No clear tendencies can be inferred, except that exhibits some kind of "plateau".

The critical angle θ c separating the two regimes is almost the same (close to 30 in Sec. 3.2.2 that the experiments performed by ONERA within the framework of PHYSICE2 have indicated that the value of ∞ depends strongly on the wall characteristics. Figure 3.6 from [1] shows a comparison for the collection efficiency β dep between experiments and numerical results from NASA and ONERA (denoted here "present SLD model"). The "full deposition" (β dep = β imp ) model is shown as well. A trend which is observed and not yet explained is that the model prediction is better on the suction side (s > 0) than on the pressure side (s < 0). This could be a guideline for considering local aerodynamic effects in the model for . In the same way, the influence of the SLD model has been studied on the final ice shape (Fig. 3.7). The conclusions are similar to those illustrated in Fig. 3.6, namely:

• The new SLD model leads to a significant improvement of the collection efficiency prediction compared to the full deposition model.

• Impingement limits (which correspond to small values of θ) are accurately considered by the SLD model.

Supplementary SLD experiments were put in place at ONERA (Sec. 3.2.2). The objective was to provide a direct measurement of the sticking efficiency for a wider range of droplet velocities, impact angles and diameters than in the Papadakis et al. database.

Extension of the Papadakis et al. database. Studies conducted by ONERA in the framawork of the project PHYSICE2

Within the framework of the DGAC PHYSICE2 convention, new experimental databases were provided:

• BASE1: a wide range for the angle of incidence θ was swept. This covered grazing impacts (θ = 10 • ) to normal impacts (θ = 90 • ). However, the velocities ranged from 25 m.s -1 to 56 m.s -1 , which was below the velocities used in the Papadakis experiments (78 m.s -1 ).

• BASE2: compared to BASE1, the droplet impact velocity was increased between 40 m.s -1 and 81 m.s -1 . The minimum observed angle was θ = 20 • . The objective was to focus on normal impacts at high velocity and to study the influence of the wall coating on the splashing phenomenon. The wall was a polished aluminium plate.

• BASE3: as in BASE2 except that the wall was covered with blotter paper similar to that used in Papadakis' experiments.

Comparisons were made with Papadakis' databases. Figure 3.8 shows the sticking efficiency (Eq. (3.4)) as a function of the incidence angle θ. The databases of Papadakis, BASE1, BASE2 and BASE3 are superimposed. In addition to the symbols representing the sticking efficiencies measured either experimentally for ONERA bases BASE1, BASE2 and BASE3, or obtained via Eq. (3.4) for the Papadakis database, the convex hulls are plotted, grouping impacts at low angles of incidence (dashes) and those at higher angles of incidence (continuous lines). Two regimes, which were already identified in [1] with the Papadakis databases (Sec. 3.2.2), are also found with all the ONERA databases (Fig.

3.8)

. First of all, for grazing impacts, the droplet sticking efficiency is an increasing function of the angle of incidence. The second regime concerns normal impacts where the sticking efficiency reaches an asymptotic value at high velocities. This asymptotic value depends on the wall characteristics (covered or not with blotter paper). Detailed conclusions are provided in [12]. In addition, PHYSICE2 was the opportunity to derive correlations for the characteristics (diameter, velocity and direction) of the secondary (re-emitted) particles. The impingement limit between bouncing and splashing regimes was studied too. These works are summarized in the PHYSICE2 synthesis report [12]. (Papadakis). Convex hulls are drawn for small angles of incidence (dashes) and for higher angles of incidence (solid lines). There are two regimes: a regime where tangential effects are predominant and a regime where normal effects are predominant. Reprinted from [12].

Direct Numerical Simulations (DNS) of the impact of a single droplet

This section summarizes the PhD works of T. Xavier I co-advised with D. Zuzio and J.L. Estivalezes. They are presented in [6]. An article on this topic is being written. Context is DNS of the droplet impact on a wall covered or not by a liquid film. DYJEAT, the code used to perform the simulations, is described in Sec. 8.2.

Figure 3.9 represents the droplet impact on a dry wall at low velocity (∼ 1m.s -1 ). DNS are compared to experiments from Yokoï et al. [44]. Qualitatively, the computational simulations are in good agreement with the experiments at any time. To achieve this, a model for the dynamical contact angle (see Sec. 7) was chosen (Fig. 3.10a), namely the one from Yokoï et al. [44] with a static angle of 90 • . A comparison between the droplet spreading diameter computed by DYJEAT and the one obtained from the experiments of Yokoï et al. [44] is presented in Fig. 3.10b.

The droplet impact on a dry wall at higher velocity (∼ 30m.s -1 ) is shown in Fig. 3.11. The qualitative agreement between DNS and experiments is rather good. However, it is worth noting that the ρ water /ρ air density ratio has been reduced to 100 (instead of 1000 in the experiments) to be able to get the corona splashing. DNS of the droplet impact on a dry wall appeared to be very challenging. Indeed, at the very first moments of the impact, lamellas are created and skate on a film of air which is so thin that it appears difficult to be meshed, even with a DNS approach. In addition, the supposed velocities in this film of air are in the order of 100 m.s -1 , so that the assumption of incompressibility may be compromised. Moreover, taking into account the dynamics of the contact line is an additional difficulty.

For all the aforementioned reasons, the droplet impact on a wall covered by a liquid film has been Reprinted from [6].

droplet impact for the Cossali et al. configuration [45].

To summarize the current work, studies have been so far focused on the normal impact of droplets. When the wall is dry, DNS struggle with capturing the triggering of the splashing phenomenon. Tricks like decreasing the density ratio are then necessary to observe droplet breakup [46]. It appears that many physical phenomena related to the droplet/wall interaction (contact line modeling) or to the wall itself (surface texture, roughness, presence of a liquid film) and which are inputs for DNS, are in fact not well controlled. Focus will now be given to the more realistic impacts on a wet wall. Preliminary computations have been initiated and tend to prove that the water film thickness plays a major role in the dynamics of splashing. Oblique impacts will also be studied. Regarding the objective of deriving a model from DNS for the sticking efficiency and for the characteristics of the re-emitted droplets, it appears to be too ambitious for the moment. Indeed, a parametric study based on the droplet diameter, velocity, angle of incidence θ is out of reach due to high CPU costs. The relevance of DNS is rather in the understanding of small-scale physical phenomena at early stages of splashing: initiation and growth of the lamella, air cushion dynamics, crown digitation . . .

Conclusions and perspectives

The wide dispersion of the results in the databases has not allowed to derive a model compatible with all databases. However, common trends with two regimes have been identified. First of all, for grazing impacts, the droplet sticking efficiency is an increasing function of the angle of incidence. After analysis of the correlation coefficients, the parameter that seems most relevant to describe sticking efficiency at low incidence angles is the angle of incidence itself. The second regime concerns normal impacts where the sticking efficiency reaches an asymptotic value. The impact velocity seems to be the most relevant parameter to describe sticking efficiency for this regime. At high velocities, the asymptotic value of the sticking efficiency also depends widely on the surface properties.

DNS of the normal impact of a SLD have demonstrated an interest for the understanding of basic mechanisms driving splashing, such as the presence of a liquid water film on the wall or the influence of pressure. However, the derivation of a macroscopic model from DNS for wall sticking efficiency is not yet possible. New approaches are being investigated, such as the impact of grazing droplets and the consideration of the solid phase in addition to the two other ones (gas and liquid) using immersed boundary methods (IBM) for instance.

Investigations on the sticking efficiency for SLDs will be continued within the framework of the European project ICE GENESIS. For this purpose, the icing wind tunnel developed at ONERA will be used. The European projects HAIC (completed in 2016) and MUSIC-haic (started in 2018) have been the frameworks of my involvement in the development of models and numerical tools for Ice Crystal Icing (ICI). The models are described in Secs. 4.2, 4.3 and 4.4 and are the subjects of 4 papers [1,2,3,4] and several conference papers or oral presentations [5,6,7,8,9,10]. These projects have been the opportunity for me to enforce relationships with some international partners (writing of common articles for instance).

Chapter 4

Ice crystal icing

Introduction

Commercial aircraft have been experiencing in-service events (probe blockage, engine rollbacks and flameouts . . . ) while flying in the vicinity of deep convective clouds since at least the early 1990s. Information gathered on over 100 weather related engine power loss events [11] has permitted the Scientific and Regulatory community to conclude that aircraft flying through areas of high Ice Water Content (IWC) are subject to a specific type of weather induced incidents. High water content is often found in deep convective clouds present in the warm tropical regions. These clouds can contain deep updraft cores that transport low-level air high into the atmosphere, during which water vapour is continually condensed as the temperature drops. In doing so, these updraft cores may produce localized regions where very high concentration of ice crystals can be encountered and where ice particles may also be found simultaneously with supercooled droplets (mixed phase icing conditions). Compared to the case of supercooled water icing conditions (Sec. 3), major modifications have to be introduced as far as the modeling of ice-crystal and mixed-phase ice accretion is concerned [5]:

• Since ice crystal accretion is supposed to mainly occur in a warm environment, it is mandatory to consider melting of the ice particles, evaporation / sublimation and heat exchange with the air flow along their trajectory.

• Ice crystals being non-spherical, the influence of their shape must be taken into account in the expression of the heat and mass exchange coefficients and also in the expression of the force exerted by the air on the particles.

• Since ice particles (or partially melted ice particles) may bounce, shatter and partially stick when they impinge on a dry or a wetted surface, an impingement model must be introduced which takes into account the influence of wall temperature, film thickness and particle liquid water mass fraction on the outcome of the impact.

• Re-emitted particles must also be accounted for since they may re-impinge downwind in the engine core.

• As shown by experimental results, erosion effects may have a strong influence and must also be considered in the impingement model.

• For ice crystal icing conditions, there are experimental evidences that liquid water may be trapped among the accreted ice particles, leading to the formation of slushy ice deposits with a porous media.

• The presence of ice particles (or partially melted ice particles) must also be taken into account in the Messinger mass and energy balance equations.

• Last but not least, in the primary stages of ice crystal accretion, heat exchange with the wall and unsteady effects must also be taken into account, which is not the case in the classical Messinger model [12].

To observe ice growth on a surface due to ice crystals, liquid water has to be present to allow ice particles to stick (at least partially) to the wall. This condition is obviously met in mixed phase clouds (i.e. a mix of solid ice crystals with liquid droplets, see hereinbelow) but also in high altitude glaciated clouds with warm environment( heated probes or engines) where initially fully frozen particles can melt. For engines the following scenario, as proposed by Mason et al. [11] has been assumed:

• The warm airflow inside the engine melts some of the ice crystals

• When droplets and partially melted ice particles impact a solid surface, a water film can be created on the surface.

• The water film captures incoming ice particles long enough for heat transfer to take place.

• Heat is extracted from the surface until the freezing point is reached, and ice begins to form.

• After this point, it is contended that further impingement of liquid and ice particles on the surface would accrete and shedding could occur.

Heated probes and engines are the areas of aircraft most prone to mixed phase and glaciated icing threat. To address the points raised above and in anticipation of regulation changes regarding mixed phase and glaciated icing conditions, the European HAIC project was launched in 2012. Its objective was to provide the necessary Acceptable Means of Compliance (numerical and test capabilities) and appropriate ice particle de-tection/awareness technologies to the European aeronautical industry in order to enhance safety when an aircraft is flying in such weather conditions. HAIC was a 4-year integrated project comprising 34 partners representing the European stakeholders of the aeronautical industry from eleven European countries and 3 partners from Australia, Canada and the United States. It was splitted into 7 subprojects. ONERA/DMPE was involved in subprojects SP5 for the development of ice crystal capabilities for ice wind tunnels (IWT) and SP6 for the development of models and numerical tools dedicated to Ice Crystal Icing (ICI). I personally took part in the activities of the SP6 where models were developed from academic experiments either proposed in the framework of HAIC or derived from literature (WP 6.3, 6.4 and 6.5). These models were implemented in the 2D icing suite IGLOO2D (Sec. 5) and assessed through a benshmark. Workpackages 6.3, 6.4 and 6.5 can be summarized as follows. More details are available in the following sections.

• WP 6.3: Ice particle trajectory model. The objective was to develop (statistical) numerical models for ice-particle trajectories which address both dynamic and thermal aspects (see Sec. 4.2). Our works were published in two papers. In [1], in collaboration with T. Hauk and E. Bonaccurso from Airbus Group Innovation, a model for ice crystal melting is proposed. In [2], in collaboration with E. Norde from the University of Twente, a comparison between an Eulerian and a Lagrangian approach for simulation of ice crystal trajectories and impact in a generic turbofan compressor is provided. The evolution of the models as well as the different validations and assessments were proposed in several conference papers or oral presentations [5,7,8,9,10].

• WP 6.4: Ice particle impingement model. The objective was to develop numerical models for computing ice particle impingement on a dry or wet surface and consider influencing parameters such as particle size, velocity and temperature, film thickness, wall temperature or presence of liquid water around the particle ice core (see Sec. 4.3). Two articles were written. A model is proposed in [4] for the glaciated regime from the experiments of T. Currie from the NRC. Another model dedicated to the mixed phase regime is proposed in [3] from a collaboration with A. Baumert and S. Bansmer from T.U. Braunschweig. In [2], realistic wall interaction models predicting rebound, shattering or deposition are applied. The different versions of our models were also proposed in several conference papers or oral presentations [5,6,9].

• WP 6.5: Ice particle accretion model. The objective was to develop and validate numerical models for computing ice accretion phenomena and predicting ice growth on a surface impacted by ice crystals (see Sec. 4.4). To do this, the classical Messinger model [12] was extended to ice crystal icing. Notably, we proposed a revisited model for runback in [4,3]. Other more or less conclusive formulations of the models were also proposed in several conference papers or oral presentations [5,6,9]. WPs 6.6 and 6.7 were dedicated to the industrial rollout of the models and tools developed in WPs 6.3, 6.4 and 6.5.

Ice particle trajectory models

Ice particle trajectory models are detailed in [1,5] and reminded in [9].

Geometrical parameters and particle motion model

As shown in Fig. 4.1, the natural and artificial ice crystals may have very complex and irregular shapes: aggregates, plates, columns . . . The maximum dimension of different ice particles varies August 22, 1999 [13]. Reprinted from [1].

between several microns and a few millimeters (see the introduction of [1] for a complete description of the measurement campaigns dedicated to the collection of data about ice particles in deep convective clouds). Moreover, fast impact of an ice crystal on a wall may also produce fragments with a complex shape. Therefore, given the wide variability in the ice crystal shapes and for modeling purpose, it is easier and more general to take into account the shape influence by introducing some global geometric descriptors:

• The particle equivalent diameter d p = 2r p = 6Vp π 1/3 • The particle sphericiy Φ = πd 2 p A • The particle crosswise sphericity Φ ⊥ = πr 2 p A proj,⊥
where A is the particle total area and A proj,⊥ the flow projected area of the particle. Using these global parameters, the drag force is given by:

F D = - 1 2 ρ a C D A proj,⊥ |w|w (4.1)
where ρ a is the air density and w = v p -v a is the relative particle velocity. It may be rewritten:

F D = - π 8 ρ a C d d 2 p |w|w (4.2)
where by definition, the normalized drag coefficient C d is defined as:

C d = C D Φ ⊥ (4.3)
Several expressions can be found in the literature for C d as functions of Φ, Φ ⊥ and Re p the Reynolds number based on the particle equivalent diameter d p . We have used the Clift and Gauvin like models of Heider and Levenspiel [14], Ganser [15] and Hölzer and Sommerfeld [16]. These models were assessed during the HAIC project, in particular for the TRL4 benchmark on trajectory models [17] where the trajectories of various non-spherical particles carried in non-uniform air flow around a cylinder were captured using high-speed videos. The agreement between the numerical and the experimental paths is globally very good even if the available experimental results did not permit to establish a hierarchy between the three semi-empirical drag models found in the literature. However, for highly non-spherical particles (with aspect ratio very different from 1), the classical Schiller and Naumann's drag model with no shape correction has been shown to be really not suitable and could lead to important errors.

Heat exchange and phase change model

As mentioned above, ice crystal icing is supposed to mainly occur in a warm environment. Therefore, the thermal aspects of an ice particle along its trajectory (melting, evaporation/sublimation, heat exchange with surrounding hot air) have to be accounted for. In the past decades, only a few studies have been carried out to improve the understanding of the melting process of individual ice particles in an airflow. Among them, few concern the melting process of individual nonspherical ice particles whose shapes are similar to the ones that can be found in the vicinity of deep convective clouds. Moreover, the particular configuration regarding the melting of ice particles in a jet engine compressor occuring at ambient temperatures (more precisely at wet bulb temperatures, see Secs. 4.3 and 4.4) higher than 0 • C has never been studied. Past works include that of Rasmussen et al. [18,19,20] focused on ice spheres or that of Matsuo and Sasyo [21,22], Fukuta et al. [23] or Mitra et al. [24] focused on the melting of snowflakes. The authors developed a theoretical heat transfer model for snowflakes made of dendritic snow crystals. The results agreed well with experimental data. However, once again, the considered shapes were far from those observed for ice crystals in convective clouds. We proposed a heat exchange and phase change model in [1] for spherical and nonspherical particles. It is as well described in [5,9]. The theoretical melting times obtained for ice particles with this model are compared to the experimental ones. A similar model was proposed in [25] where the capillary flow of the melted liquid around the solid particle as well as the effects of the interface curvature are taken into account. The initial shape of the particle in [25] is approximated by a spheroid and the interface of the particle is tracked by a level-set method. Apparently more sophisticated, this model [25] however gives results similar to our model [1].

Heat exchange model

In absence of phase change phenomena (sublimation, evaporation, condensation, melting), the general form of the heat equation reads [5,9]:

m p c p dT p dt = πd p N u Φ k a (T a -T p ) (4.4)
where m p is the particle mass, c p the specific heat capacity of the particle, T p the particle temperature (supposed to be almost constant), T a the air temperature, N u the Nusselt number and k a the air heat conductivity. A model for N u extended to nonspherical particles was developed. In case of spherical particle, the most commonly used model for the Nusselt number is the Frossling correlation [26] which reads: N u = 2 + 0.55 Re p P r 1/3 (4.5)

In [5,1], from a Reynolds-type analogy, we proposed the following generalization of Eq. (4.5):

N u = 2 √ Φ + 0.55 Re p P r 1/3 Φ 1/4 (4.6)
Validations of this model were performed in [5] from available results in the literature: the ones from Richter et al. for prolate spheroids and cuboids [27] and the ones from Comer for oblate spheroids [28] and P r = 0.7.

Phase change model

Modeling of the melting process was presented in [1,5]. It is assumed to be divided into three successive phases (Fig. 4.2):

• First phase (Fig. 4.2, left): the ice particle is solid and is heated from its initial temperature to the melting temperature equal to 0 • C. A mass transfer occurs because of sublimation and deposition. Ice particle temperature is supposed to be uniform. Its density is constant as well as its sphericity Φ.

• Second phase (Fig. 4.2, middle): the particle starts melting at constant temperature T p = 0 • C until all solid ice becomes liquid water. A mass transfer into and from the environment occurs due to evaporation and condensation. The ice core is completely surrounded by a liquid film of uniform thickness. Particle sphericity progressively increases up to 1 according to the following law:

Φ = (1 -η m ) • Φ 0 + η m (4.7)
where Φ 0 is the solid ice crystal sphericiy before it begins to melt and η m is the ice crystal melting ratio (i.e. the liquid water mass fraction present in the melted ice particle). The liquid layer is supposed to be at rest. Shedding of the melt water does not occur. Density of the melted ice particle is given by:

1 ρ = (1 -η m ) • 1 ρ i + η m • 1 ρ l (4.8)
with ρ i and ρ l respectively the densities of the ice core and the liquid water. • Third phase (Fig. 4.2, right): the particle is a pure liquid droplet and is heated from the melting temperature 0 • C. A mass transfer occurs because of evaporation and condensation. Droplet temperature is supposed to be uniform. Its density is constant as well as its sphericity Φ which is equal to 1.

The equations governing each phase are presented in [1]. For the mass transfer equations, the particle Sherwood number Sh is obtained from an analogy between heat and mass transfer:

Sh = 2 √ Φ + 0.55 Re p Sc 1/3 Φ 1/4 (4.9)
where Sc is the Schmidt number.

Comparisons between theoretical and experimental results

Experimental works dedicated to the study of the melting process of ice particles were driven by T.

Hauk from Airbus Group Innovation. The experimental apparatus consists of an acoustic levitator, aimed to suspend individual ice particles; an optical and recording system to magnify and record the melting particle; and a flow control system to generate a defined airflow. See [1] for a complete description. The experimental results were compared to the ones obtained from the models described previously. In total, 222 melting processes were recorded. The airflow parameters were varied to reproduce 13 different conditions at which the experimental (and theoretical) melting times were obtained. At each condition, a small number of spherical ice particles were melted, followed by several melting processes of nonspherical particles. The different test conditions are listed in Tab. 1 from [1]. The initial ice particle temperature, the ambient pressure and the flow velocity varied little from condition to condition. On the other hand, for each condition, great differences were observed for the flow temperature and the relative humidity, which allowed a wide wet bulb temperature range to be scanned (from 6 to 20 • C). The maximum dimension of the nonspherical ice particles placed in the acoustic levitator varied between 730 and 3490 µm.

For the validation of the melting model, the almost-spherical particles of each condition, which exhibit sphericity values close to 1, are considered (Fig. 4.3a). In general, there is an excellent agreement between the theoretical and experimental melting times for the considered ice spheres. In only 1 out of 65 melting processes, the experimental melting time is significantly larger than the theoretical melting time. For this particular case, the error bars of both times do not overlap. The mean value of the relative differences is 8.1%.

In Fig. 4.3b, the experimental melting times of 157 nonspherical ice particles are shown and compared with the theoretical melting times of ice spheres obtained with the assumption Φ = 1. In Fig. 4.3c, the corresponding histogram of the relative differences is shown. Theoretical melting times tend to overestimate experimental melting times. This overestimation is significant in 31 out of 157 cases. The mean value of the relative differences is 19.9%.

In Fig. 4.4a, the experimental and theoretical melting times calculated with the assumption Φ = Φ ⊥ are shown (nonspherical ice particles). In Fig. 4.4b, the corresponding histogram of the relative differences is shown. Conditions 1 (relative humidity of 4 %) and 4 (relative humidity of 74 %) are shown separately in Fig. 4.4a to outline the effect of evaporative cooling. This confirms the fact that evaporative cooling plays a major role in condition 1 and increases the melting times significantly because of the low relative humidity of the flow. A very good agreement between the theoretical and experimental melting times is reached. Only 5 out of 157 nonspherical particles show significantly larger or smaller theoretical melting times than measured. The mean value of the relative differences is reduced to 9.7 % (compared to 19.9 % when Φ = 1 is assumed in the models).

These experiments have validated the heat exchange and phase change models. For non-spherical ice particles, it was demonstrated the importance of taking into account the shape of the particle (via the sphericity). The first order dependency of the models to particle sphericity explains all the difficulty to evaluate accurately the melting time for an ice crystal. Indeed the sphericity of an ice particle can, at the best, be approximated by the experimental measurement of the orthogonal sphericity. Therefore, the uncertainty on the measurement of Φ has a direct impact on the particle melting time. Another key parameter is the evaporative cooling effect (or the relative humidity of the surrounding air). For instance, this parameter proved to be of paramount importance in assessing characteristic melting times for ice particles in the NRC RATFac facility. Indeed, the configuration where cold (≈ -15 • C) and dry air loaded with pure ice particles is mixed with warm moist air to control the humidity level and the melting ratio of the ice particles is met.

Ice particle impingement models

We propose here a description of the models developed in the framework of HAIC and dedicated to ice crystal impingement as well as the experimental data used to derive these models. For a detailed description, see [29,30,5,9]. Due to the lack of experimental data, it has been very difficult to develop a comprehensive model which accounts for ice particle / solid surface and ice particle / liquid film interactions. The laboratory experiments were performed by T. Hauk [31,32,33] in order to fill existing gaps and to provide experimental data for model calibration and validation. However, note that a limited range for the impact velocity (V imp ≤ 74 m/s) was observed in the experiments from T. Hauk, which is far from the real operating range (V imp up to 300 m/s). The models presented in [9] are a refined version of the ones proposed in [5]. They are based on Hauk's experimental results, literature data [34,35,36,37,38,39] and several simplifying assumptions. The models for particle impingement can be divided in two subcategories:

• Models for the sticking efficiency. When they hit a wall, the particles may stick or partially stick to it. Models are needed to evaluate the mass fraction of the ice particle which remains on the wall after the impact.

• Models for the secondary particles. After a wall impact, depending on the impinging particle characteristics and wall properties, it is assumed that ice crystal impingement may lead to three possible outcomes [9]:

-Sticking regime: the particle totally sticks to the wall ;

-Bouncing regime: the particle bounce off the wall with partial deposit. Particle velocity is modified but its size and shape do not change ;

-Shattering regime: the particle shatters into several small fragments. Some fragments are re-emitted into the air flow and the others stick to the wall during the impact (due to the presence of a liquid film on the wall or the presence of liquid water inside the impinging ice particle itself).

Models for the re-emitted particle diameter, melt ratio, ice core density, sphericity and velocity (or, equivalently for the restitution coefficient of the impinging kinetic energy) have to be proposed.

I have been personally involved in the development of the models for the sticking efficiency [5,9,3,4]. On the other hand, regarding the models for the secondary particles, this has not been the case. I just integrated the models proposed in the framework of HAIC into the simulation tool IGLOO2D (Sec. 5) . Indeed, the main concern for secondary particles is re-impingement which appeared to be of second order for the academic geometric configurations used to validate the models [3,4]. Obviously, this is not the case for confined geometric configurations such as those encountered in aircraft engines or in high lift configurations for a wing. We (at ONERA) initially decided to focus on the modeling of the sticking efficiency, the erosion rate and the liquid water mass fraction inside the slushy deposit (Sec. 4.4).

Regarding impingement models, heat exchanges between the wall and the ice particle during the impact have been supposed to play a negligible role. This assumption is not restrictive as far as accretion on unheated surface is concerned (except during the short transient phase during which the wall temperature is decreasing). This is the hypothesis assumed in Sec. 4.3.1. For a heated surface, there exist experimental evidences showing that particle may partially melt upon impact and that a small amount of water may stick to the wall, promoting further ice crystal accretion. Experiments were dedicated to the impingement of ice particles on heated surfaces and are briefly described in Sec. 4.3.2. However, contrary to unheated surfaces, no universal models have been derived. They will be considered in forthcoming works (MUSIC-haic).

Impact on unheated walls

Models for the secondary particles

As explained previously, I have not been directly involved in the development of the models for the secondary (re-emitted) particles. I have just participated in the discussions within the framework of HAIC and integrated them in IGLOO2D. However, for completeness, I briefly summarize them here. For the simplified configuration of a pure ice crystal impacting on a dry wall, the onset between the bouncing and the shattering regimes strongly depends on the particle normal kinetic energy before the impact [39]. To account for this experimental fact, the L number introduced by Vidaure and Hallet [39] for ice crystals and by Guégan et al [37] for hail particles, is used:

L = πρpd 3 p (vpn) 2 12 πe σ d 2 p = 1 12
ρ p d p (v pn ) 2 e σ (4.10)
L corresponds to the ratio of the normal kinetic energy to the surface energy (assuming a spherical shape for the particle). It can be compared to the Weber number for SLD (Sec. 3.2) where e σ is the energy per unit area related to the creation of cracks inside the particle. e σ is the counterpart of the surface tension coefficient σ for the droplets. An expression for e σ is proposed in [5]. The transitions between bouncing (low values for L) and fragmentation (large values for L) are modeled very simply [9]:

• L ≤ L c1 : elastic rebound without particle fracturing and kinetic energy loss (restitution coefficient equal to 1).

• L c1 ≤ L ≤ L c2 : inelastic rebound with plastic deformation, particle internal fracturing and loss of kinetic energy (restitution coefficient lower than 1).

• L c2 ≤ L: highly inelastic impact with major particle breakup, important loss of kinetic energy (restitution coefficient much lower than 1) and secondary particle re-emission.

The values for the parameters L c1 and L c2 were derived from the experiments of T. Hauk [30].

Regarding the bouncing regime, if ξ B t and ξ B n respectively correspond to the tangential and normal restitution coefficients of the kinetic energy, and in the absence of experimental data for ξ B t , it is assumed that ξ B t = 1 which is in agreement with the experimental observations for hail particles [36]. For ξ B n , the results obtained by T. Hauk in the framework of HAIC [30] are used. Figure 4.5a shows a comparison between the model for ξ B n , the experimental results from T. Hauk and the results from Pan & Render [35] for hailstones.

Regarding the shattering (fragmentation) regime, three coefficients can be defined [9]. ξ F t is the tangential restitution coefficient of the kinetic energy, ξ F nt and ξ F nn are normal restitution coefficients. ξ F nt is the fraction of the normal momentum which is transferred into tangential momentum due to the impact. Once again, like for ξ B t , ξ F t = 1 is assumed. The model for ξ F nn is chosen similar to ξ B n . The model for ξ F nt is deduced from the experiments of Guégan et al [37] which are devoted to the impact of hail particles at high velocity. Figure 4.5b shows a comparison between the model for ξ F nt and the experiments from T. Hauk. Figure 4.5c shows a comparison between the model for the secondary (re-emitted) particle diameter [9], the experimental results from T. Hauk (AGI) and from Pan and Render [35]. The agreement is reasonably good over a wide range of values for the impact parameter L.

For information, TU Darmstadt has developed an impact model for ice crystals in the framework of HAIC [30]. The model of the particle collision with a dry rigid substrate and breakup is capable of predicting the shattering of the ice crystal depending on its size and velocity. It is based on a physical consideration of the occurring stresses in the particle [40,31,33].

The impact of a rigid sphere on a free liquid water surface is considered in [41,42,43]. The agreement between numerical prediction and experimental data is very good. However, the ratio between the thickness of the liquid water film and the particle diameter is much higher than that found in realistic configurations encountered in flight. T. Hauk [31] presented as well an experimental method dedicated to the observation of spherical ice particle impacts onto thin liquid water films in a subfreezing environment. He finds that for relatively thin dimensionless water films (0.05 ≤ h/D 0 ≤ 0.25) the maximum no fragmentation velocity is approximatively twice as high as in the case of a dry wall.

Models for the sticking efficiency

One of the main challenges has been to estimate the deposited mass rates of liquid water ṁL dep (coming from both impinging liquid droplets and melted ice crystals) and of solid ice ṁS dep (coming from impinging ice crystals). They are defined as:

ṁL dep = ṁimp,d + s η m ṁimp,c ṁS dep = s (1 -η m ) ṁimp,c (4.11)
where s is the ice crystal sticking efficiency and ṁimp,d (resp. ṁimp,c ) the impinging droplet (resp. ice crystal) mass rates. s correlates with the difference of impinging and rebounding particles. η m is the ice crystal melting ratio. To set the notations, we can write: where ṁL imp,c and ṁS imp,c are respectively the ice crystal liquid and solid impinging mass rates. With no liquid water, the ice particles do not adhere to cold airframe surfaces and bounce off. Regarding the presence of liquid water which is necessary for the ice crystals to stick, two origins are possible. On the one hand, the liquid water may come from supercooled water droplets from mixed phase clouds at atmospheric temperatures above -40 • C. These atmospheric conditions are referred to "mixed phase conditions" and are encountered at temperatures below 0 • C. On the other hand, the liquid water may come from the melted part of the ice crystals themselves. This regime is referred to "glaciated conditions" and is often associated to engine conditions characterized by a wet bulb temperature above the freezing point. Both glaciated and mixed phase conditions occur in convective clouds and have been present during engine power-loss and damage events (Mason et al. [11]). Generally speaking, the amount of liquid that is involved in ice crystal icing is characterized by the liquid water ratio LW R, which is defined as the ratio of liquid to total water content:

ṁL imp,c = η m ṁimp,c ṁS imp,c = (1 -η m ) ṁimp,c (4.12)
LW R = η m IW C + LW C T W C (4.13)
where LW C is the liquid water content, IW C the ice water content, η m the ice particle melting ratio and

T W C = LW C + IW C.
The modeling of the general scenario, i.e. the mixing between melted ice crystals and supercooled water droplets, has turned out to be too far ambitious. That is why we have been focusing on the modeling of two simplified configurations:

• The glaciated conditions with no supercooled liquid droplets. In Eq. (4.13), LW C = 0, LW R = η m and the liquid water comes only from the melted part of the ice crystals themselves. This configuration was published in [4] and was largely based on T. Currie's experiments.

• The mixed phase conditions where the liquid water comes from supercooled water droplets only. The ice crystals are pure solid cores. η m = 0 and LW R = LW C/T W C. This configuration was published in [3] and was largely based on the experiments from A. Baumert (T.U. Braunschweig).

Regarding glaciated conditions, we proposed in [4] the following model for s (Eq. (4.11)):

s = F (η m ) (4.14)
It should be noted that this model is intentionally simplified in the sense that it only takes into account the most relevant parameter, namely η m . However, it can be expected that a more comprehensive model would be a function of other parameters such as the angle of impact or the magnitude of the impact velocity. The function F has to meet the following conditions:

1. In the absence of liquid water, pure ice crystals (η m = 0) bounce off the wall ( s = 0). Therefore, F (0) = 0.

2. For the supercooled water icing regime (η m = 1), all the droplets stick to the wall ( s = 1). Therefore, F (1) = 1.

3. In the experiments from Currie et al. [44], it was shown that for a crowned cylinder, the sticking efficiency near the stagnation point is a function of the particle melting ratio η m . More specifically, for low η m , s s is a linear function of η m (see hereafter). Therefore, F (η m ) ≈ K # • η m at low η m (with K # an adjustable parameter to be calibrated).

4. The function F is expected to be a smooth increasing function which tends to 1 when η m → 1.

Given all these conditions, the following polynomial expression is chosen for F [4]:

F (η m ) = K # -2 η 3 m + 3 -2K # η 2 m + K # η m (4.15)
In the experiments of Currie [44], the apparent local sticking efficiency actually results from the competition between two independent phenomena: the accretion of new incoming ice crystals and the erosion of the accreted ice layer by impinging particles. At numerical level, each of these phenomena can be described separately. This is not the case for the experimental observations where only the apparent net sticking efficiency is measurable, which combines both the partial sticking of the ice particles and the erosion effects. However, the erosion rate depending mainly on the tangential impact velocity, it is minimum near the stagnation point where the velocity is mainly normal to the wall (see Sec. Reprinted from [4].

Regarding mixed phase conditions, we proposed in [3] the following model for s (Eq. (4.11)):

s = F (LW R stick ) (4.17)
where LW R stick is the ratio of liquid to total water content at the wall necessary for the ice crystals to stick and F is the function defined in Eq. (4.15). Just above, for glaciated conditions, we have supposed that LW R glaciated stick = η m . For the mixed phase regime, many attempts have been made to model LW R mixed stick . In a first version of our model [5], LW R mixed stick was supposed to be a function of the wall water film thickness h w . However, h w is very difficult to be accurately measured especially when it is in the order of a few tens of microns. A wrong estimation of h w has turned out to have a strong influence on the final ice shape. To avoid these pitfalls, we proposed in [9] a model for LW R mixed stick defined as a function of the ratio LW C/T W C. Finally, it appeared that the accretion surface total liquid mass fraction f l (Eq. (3.2)) was more relevant than the classical LW C/T W C proposed in [9]. Indeed, for pure ice crystals where LW C/T W C = 0, the only available liquid water which allows ice crystals to stick to the wall comes from the liquid film on the wall which is itself defined by f l > 0. Finally, we proposed for LW R mixed stick [3]:

LW R mixed stick = K d • f l (4.18)
where K d is an adjustable parameter calibrated in [3].

Note that today, a general model for s in the case of melted ice crystals mixed with supercooled water droplets does not exist. One of the reasons is that we have not yet really understood how high the liquid water encountered in the ice crystal itself (via η m ), in the supercooled water droplets (via LW C) or at the wall (via f l ) plays in the sticking efficiency coefficient s . Clarification will have to be provided in MUSIC-haic (see Sec. 4.6).

Impact on heated walls

Ice particle impingement on an heated wall was studied in the framework of HAIC. For instance, in TsAGI experiments [30], the interaction of ice crystals with a heated airfoil leading edge was investigated. Experiments were performed for both a dry heated wall and a heated wall with a liquid film. The liquid film was created by a swirl water injector located upstream the model. The crystals rebound/sticking threshold was investigated by varying the wall temperature or by varying the crystals diameter while the other parameters were kept constant.

The models developed at the ONERA in the framework of HAIC did not take into account the effects of heated walls on ice particle impingements. This will be the subject of specific studies within the framework of MUSIC-haic (see Sec. 4.6).

Ice particle accretion models

Ice crystal icing is a combination of several physical phenomena. Firstly, ice crystals of different shapes, sizes and velocities impact onto a solid substrate (Sec. 4.2). Some of them may stick to the substrate provided that there is some liquid water. The ice particle deposited mass rate is determined by the sticking efficiency coefficient s (Sec. 4.3). Secondly, the ice layer over the substrate is a porous media. Therefore, a part of the liquid water remains trapped and accumulates in the pores of the ice layer. This leads to a new formulation of the runback model for the liquid water. Finally, a part of the accreted ice and water is ejected by the impacting ice crystals. This phenomenon is described by the erosion model. In this section, the extended Messinger balance for ice particles is presented as well as the new formulation of the runback model (Sec. 4.4.1). Finally, the erosion model is described (Sec. 4.4.2). These models were presented in [5,9,4,3]. This section is largely inspired by [4].

The extended Messinger balance

The classical Messinger's mass and energy balance equations [12] have been adapted to take into account the presence of ice crystals among the impinging particles. The corresponding models are accurately described in [45] and [5]. To introduce the notations, the general form of the mass balance equation is given:

ṁS acc + ṁL M ess = ṁrbi + ṁdep -ṁevs -ṁer ṁL M ess = ṁL acc + ṁrbo (4.19)
where ṁrbi is the incoming runback liquid water mass rate, ṁdep is the particle deposited mass rate, ṁevs is the evaporated/sublimated mass rate, ṁL M ess is the net liquid water mass rate (i.e. the mass rate corresponding to the liquid water which does not evaporate, sublimate or freeze during the corresponding time step), ṁS acc is the accreted ice mass rate and ṁer is the erosion mass rate. The mass rate ṁL M ess is splitted between ṁL acc and ṁrbo which respectively stand for the mass rate of liquid water which remains trapped inside the porous ice layer and the mass rate of running back liquid water. All the terms are summarized in Fig. 4 acc are respectively the mass rates for impinging ice particles, deposited ice particles, evaporation/sublimation, erosion, incoming runback liquid water, outgoing runback liquid water, liquid water trapped in the pores of the accreted ice layer and solid ice accretion. Reprinted from [4].

impinging mass rates and is described in Eq. (4.11).

In the classical Messinger model [12], the hypothesis that the liquid water may not locally accumulate and is carried downstream by the air flow ( ṁL acc = 0 and ṁrbo = ṁL M ess ) may be justified in the case of supercooled water icing where the frozen droplets form a solid non porous ice layer. But, for ice crystal icing conditions, there are experimental evidences [46,44] that liquid water may be trapped among the accreted ice particles, leading to the formation of slushy ice deposits. That is why the following runback model, i.e. the split of the water mass flow rate ṁL M ess between ṁL acc and ṁrbo , is proposed:

ṁL acc = min[H (η m ) • G (Ψ) • ṁS acc ; ṁL M ess ] ṁrbo = ṁL M ess -ṁL acc (4.20)
where G is an increasing function of the porosity of the accreted ice Ψ, which means that the more porous the ice, the more liquid water can accumulate. See [4] for more details. The function H is defined so that the classical Messinger model for supercooled water icing [12] actually appears like a degenerate configuration of the model with η m = 1 (see [4]).

Erosion model

The erosion mass rate ṁer (Eq. (4.19)) is defined as [4]:

ṁer = min ṁS acc + ṁL M ess , min (1, er ) • ṁimp,c (4.21)
The following empirical model for the erosion efficiency er is proposed [4]:

er = E V t imp,c V 0 2 1 • y l0 y l0 -min (y l , y l0 ) 2 • 1 + (l 0 κ) 2 3 (4.22)
where V t imp,c denotes the mean tangential velocity of the impinging ice crystals and y l is the wall liquid mass fraction defined by:

y l = ṁL acc ṁL acc + ṁS acc (4.23)
y l = 0 stands for a pure ice layer with no water trapped inside and which sticks firmly to the wall. y l = 1 is an asymptotic value with only liquid water. κ is the local curvature of the ice layer surface. E, V 0 , y l0 and l 0 are empirical constants to be calibrated (see the calibration section A3 in [4]). Term 1 in Eq. (4.22) is justified from experimental observations with conical ice shapes [46,44] and from the erosion effects near the impingement limits [START_REF] Currie | Fundamental study of mixed-phase icing with application to ice crystal accretion in aircraft jet engines[END_REF]. Erosion rate seems to mainly depend on the tangential component of the ice crystal impinging velocity. Term 2 defines a critical value for y l . For large wall liquid mass fractions (y l y l0 ), the slushy ice deposit at the wall contains too much liquid water and cannot stick to the wall due to erosion and continuous shedding phenomena. Term 3 accounts for the smoothing effect of erosion. The higher the local curvature, the higher the erosion rate.

Regarding the accretion models developed by the other partners in the framework of HAIC, D. Kintea from T.U. Darmstadt proposed a theoretical model for the porous ice/water layer which builds up on hot surfaces in aircraft engines or on heated probes [41,[START_REF] Kintea | Transport processes in a wet granular ice layer: Model for ice accretion and shedding[END_REF]. Numerical simulations where the individual ice particles and liquid droplets are resolved with a Volum Of Fluid (VOF) numerical method were performed too [41,[START_REF] Kintea | Numerical investigation of ice particle accretion on heated surfaces with application to aircraft engines[END_REF]. However, dynamics of particle impingement is not considered. The model accounts for phase transitions between solid, liquid and gaseous states as well as heat fluxes inside the substrate on which ice accretes and inside the porous ice/water layer. We can also mention TSAGi where a model for the hydro-thermodynamic behavior of the liquid water film streaming along the surface of a body was proposed [START_REF] Villedieu | Deliverable D65.3 -Final report on the two-layer model[END_REF].

Model validation: comparison with experimental results

The models we have developed for ICI, in particular those for the liquid water runback inside the porous media and for erosion, suffer from a high degree of empiricism. As the physics is extremely complicated and not completely understood, the proposed correlative models have been based on macroscopic experimental observations. However, and this is the point of the models we have proposed, they are able to reproduce some tendencies observed experimentally. Of course, none of these phenomena have been hard embedded in the models, but have been direct consequences observed by the application of the latter. This section is dedicated to model capabilities and assessment. Comparisons with experimental results focused on the whole icing process (i.e. from particle impingement to final ice shape) are proposed.

Regarding glaciated conditions, the models were validated from the NRC ice crystal accretion experiments in their small altitude wind tunnel [46,44]. The ice crystal melting ratio η m was adjustable by adjusting the wet bulb temperature T wb . See Sec. 4 from [4] for a complete description of the test cases used to validate the models. First of all, the numerical simulations allowed the calculation of stationary ice thicknesses whose conical shape is in agreement with that obtained experimentally (see Fig. 4.8 for the experimental shapes and Fig. 4.9b for the computational ones). Secondly, the existence of optimal icing conditions as a function of the ice crystal melting ratio (the so-called "plateau" effect) is well represented (Fig. 4.9b). This plateau (see Fig. boundary where icing severity decreases strongly. At low melting ratios (left limit) the ice particles do not contain enough liquid water to stick to the model wall. On the other hand and for large melting ratios (right limit), the amount of liquid water is larger so that the ice particles do not stick and are washed off the surface [44]. Moreover, erosion effects are strong enough to prevent ice accretion. Thirdly, the fraction of liquid water in the slushy ice layer y l predicted by the runback model is consistent with the one measured experimentally (Fig. 4.10a). Inside the ice shape as a whole, y l increases with the melting rate η m , which is consistent with an increase in the erosion rate with the melting ratio. . For all the melting ratios, the wall liquid mass fraction is greater than the melting ratio of the impinging ice crystals (y l > η m ). This tendency, which is well taken into account by the model, means that there is more liquid water in the slushy ice layer at the wall than in the partially melted ice crystals impinging the body. However, the model is not able to predict the plateau for y l obtained experimentally for η m > 0.2 (Fig. 4.10b). This is one of the weaknesses of the model. Finally, the effects of the Mach number on the ice shape is well taken into account by the model ( [4], Sec. 6.4). The ice thickness at the stagnation point is well predicted by the model when the Mach number is raised from 0.25 to 0.4. This validates the term

V t imp,c V 0 2
of the erosion model (Eq. (4.22)).

Regarding mixed phase conditions, we derived a model [3] from the experiments performed by A. Baumert from T.U. Braunschweig. During his PhD, A. Baumert completed a 6-month internship at ONERA under my supervision to use IGLOO2D. The experimental results have thus revealed shortcomings in the models from IGLOO2D that will have to be solved in future works. The first point that can be raised from Fig. 4.11 is that the sharp wedge ice shape is not accurately captured by the model. Indeed in [3] the flow velocity is U ∞ = 40 m/s while the erosion model (Eq. (4.22)) was calibrated from the NRC experimental data at U ∞ = 85 m/s with a different particle size distribution [4]. The influence of the liquid water ratio LW R (Eq. (4.13)), total water content T W C and temperature T ∞ was studied and is presented in Fig. 4.11. Regarding the influence of LW R, the numerical simulations are on the whole able to account for the global tendency (Fig. 4.11d). Figure 4.11e indicates differences in steady state accretion angle ϕ 0 depending on the melting ratio. In case of T ∞ = -5 • C the lowest steady state angle appears to be about 80 • for LW R = 0.12, indicated by the dashed line asymptote. In contrast, a rather obtuse angle of 130 • appears for LW R = 0.5. Based on accurate observations of side view video sequences [3], it is supposed that shear stresses homogenise the accretion surface by dragging the ice/water mixture in downstream direction. This kind of surface flow seems to be most dominant at high melting ratios, when significant liquid is available on the accretion surface. Thus, it is hypothesised that the obtuse angle of 130 • at LW R = 0.5 is significantly affected by a shear stress-driven transport of ice-water mixture towards the base of the accretion wedge. Such transport avoids the formation of a sharp accretion angle and reduces the accretion efficiency by conveying lose ice particles in downstream direction. This aspect is not respected in the current version of IGLOO2D. At rather dry conditions of T ∞ = -5 • C and LW R = 0.12, shear stress-driven transport is less severe and a sharp accretion angle develops (Fig. 4.11e). Figures 4.11a and 4.11b show the influence of the temperature T ∞ on the accretion shapes for two different melting ratios LW R = 0.12 and LW R = 0.28. The conclusions are similar for both LW R. IGLOO2D strongly overestimates accretion growth at the lowest temperature -15 • C because water freezing and particle sticking are not balanced properly in IGLOO2D [3]. Indeed, first, enhanced convective cooling at lower temperatures promotes a greater mass of liquid water to freeze on the accretion surface. Secondly, the temperature has an indirect influence on s , the ice crystal sticking efficiency. Enhanced freezing at lower temperatures allows less liquid water to accumulate inside the accretion layer, resulting in lower f l and less ice particle deposition (Eqs. (4.17) and (4.18)). The inaccurate ice shapes computed by IGLOO2D at low temperatures have been the main weaknesses pointed out by the experimental studies presented in [3].

As mentioned above, our models have been calibrated with the databases from T. Currie for the glaciated regime [4] and from TUBS for the mixed phase regime [3]. However, regarding mixed phase experiments, other databases have been used for validation [9]: , crystal mean mass diameter not well known (between 100 and 1000 µm).

•
The results were presented in [9]. In the framework of HAIC, the models and dedicated tools have been integrated into industrial environment, including basic validation [START_REF] Cliquet | Deliverable D66.1 -Final report on model integration in industrial environment incuding basic validation[END_REF]. Representative test cases for aircrafts, probes and engines have been included [START_REF] Cliquet | Deliverable D67.1 -Capability assessment: test case definition[END_REF].

Conclusions and perspectives

In the light of the Sec. 4.5, the models developed in the framework of HAIC are able to reproduce some tendencies observed experimentally. However, much work remains to be done to improve them and reduce empiricism. We propose now to outline some shortcomings of the proposed models that should be addressed in the forthcoming EU project MUSIC-haic. The present model is able to capture the so-called plateau effect at least for the left part of the plateau (low melting ratios). This is less the case for high melting ratios where the decrease of ice accretion is not as drastic as expected. There are two possible explanations for these discrepancies. Firstly, for high melting ratios, the adhesion of the slushy ice layer to the wall is reduced due to the high liquid water content of the ice. Thus, shedding phenomena may appear due to the aerodynamic forces. The current version of the model does not take into account the shedding process. This process should be incorporated in future versions. Secondly, for high melting ratios, ice particles do not stick and are washed off the surface since too much water is present [44]. The model for sticking efficiency (Eqs. (4.14) and (4.15)) has been calibrated for small melting ratios (see Sec. 4.3.1) and is an increasing function of η m (Eq. (4.15)). Therefore, future works should be dedicated to extend the sticking efficiency model to high melting ratios. Moreover, it was shown that the competition between water freezing and particle sticking has not been fully understood, especially at low temperatures [3].

Regarding runback, the model is not able to predict the plateau for y l when η m > 0.2. Instead, the model predicts a linear increase of y l with the melting ratio of the ice crystals. This weakness in the model is due to a poor estimation of the porosity of the slushy ice layer at the wall.

Regarding the erosion model, ongoing work is on rewriting the erosion mass rate model in terms of dimensionless numbers (using impact angle, particle kinetic energy and ice layer cohesion energy). In particular, future versions of the models for erosion and sticking efficiency will have to include the influence of the ice particle diameter since it has a significant effect on accretion [START_REF] Knezevici | Particle size effects on ice crystal accretion[END_REF][START_REF] Knezevici | Particle size effects on ice crystal accretion-part ii[END_REF]. The erosion model should take into account the shear stresses which homogenise the accretion surface by dragging the ice/water mixture in downstream direction [3]. Moreover, the normal component of the impinging velocity should be added to the erosion model (as appropriate for a brittle material) as well as to the sticking efficiency (in comparison to SLD, Sec. 3.2).

The models we proposed for ice crystal impingement have been designed for a clean unheated wall, with possible basic consideration of the water liquid film on the wall. Work is still required for heated walls covered or not by a liquid film. Impacts on walls covered by a granular ice layer is challenging too.

From all this, it was concluded that, at the end of the HAIC project, several important gaps were not solved, namely:

• Certain microphysics models are not sufficiently mature and supplementary dedicated academic experiments are needed.

• With some exceptions, the models developed in HAIC have only been integrated into 2D codes.

• To date, code validation has been limited to comparisons with academic experimental results using simple 2D geometries

The EU MUSIC-haic project, which started in 2018, is supposed to fill these gaps and incorporate in industrial 3D multi-disciplinary tools a validated ICI capability that can be used for both design and certification of new engines, probes and aircraft. I have been appointed as MUSIC-haic WP2 leader "Model development, elementary validation and down-selection". I hope that the future work and results will be commensurate with those already achieved for HAIC and that the models initiated in HAIC will reach a degree of maturity necessary to meet the objectives set out in MUSIC-haic. The first objective I had after my transfer from DSNA to DMAE has been to develop a new 2D icing suite called IGLOO2D. Its description is detailed in this chapter and was the subject of a AIAA conference paper [1].

Introduction

Several two-dimensional icing codes exist in the world: LEWICE [2,3,4] from NASA, TRAJICE2 [5,6] from DERA or ONICE2D [7] from ONERA. A comprehensive review of the models and a comparison of these codes prediction capabilities may be found in [8]. These codes all are accurate to predict ice shapes in the framework of the Appendix C [9]. From the ONERA point of view, ONICE2D was not adapted to deal with the new requirements of the Appendix C extension for SLDs (Sec. 3.2) or ice crystals (Sec. 4). A modern icing suite should be adapted to the new icing environments like engines or probes with complex geometrical configurations. Complex viscid flows at high velocities should be taken into account. Models dedicated to the complex wall/particle interactions encountered in the SLD (Sec. 3.2) or ice crystal (Sec. 4.3) regimes should be added. The accretion solver should be coupled with a film solver to improve the modeling of runback. Moreover, models for Ice Protection Systems (IPS, Sec. 6) should be added to simulate both anti-and de-icing modes (including models for ice shedding). An unsteady approach with a regular update of the ice thickness and growth direction is mandatory. Industrial partners need a modular, robust and fast icing suite. As far as these last points are concerned, three-dimensional icing codes are not yet mature enough. This explains why the two-dimensional icing codes have a place in the industrial world. For all these reasons, a new icing suite IGLOO2D has been developed at ONERA. It has been mainly funded by two French projects GENOME and CORICE and is supposed to address all the aforesaid points. In a first section (Sec. 67 5.2), an overview of IGLOO2D, mainly based on [1], is proposed. In the second section (Sec. 5.3), some test cases for the validation of IGLOO2D are presented. The framework here is the Appendix C [9]. Some references to the extension of IGLOO2D to App. O and D (respectively dedicated to SLD and ice crystals/mixed phase) are made. The reader is then referred to Sec. 3.2 and 4 respectively.

Overview of IGLOO2D

IGLOO2D is mainly composed of different stand-alone components (Fig. 5.1): The grid generation tool MESH2D is dedicated to the building of both structured (Fig. 5.2a), unstructured (Fig. 5.2b) and hybrid (Fig. 5.2c) grids. It is not limited to single profiles and grids around multi-element airfoils (Fig. 5.2d) are possible. Structured grids are computed from a home made software. Unstructured grids are computed with the GMSH [10] software. See [1] for a complete description of the mesh generation in IGLOO2D.

The solver AERO2D is supposed to compute the aerodynamic field (such as air pressure, velocity and temperature, recovery temperature, heat transfer coefficient or local steam mass fraction) necessary to the computation of particle collection efficiency and ice accretion. Firstly, the inviscid external flow is computed. Then the boundary layer equations are solved (one-way coupling). Two different approaches are proposed in IGLOO2D. The first one is based on a direct resolution of the Prandtl equations of the boundary layer [11]. The software CLICET [12] developed by the ONERA is then the suitable tool. The friction coefficient C f and the heat transfer coefficient h t are derived directly from respectively the velocity and temperature profiles inside the boundary layer. In the second approach, C f and h t derive from the integrated thicknesses of the boundary layer such as the displacement and the momentum thicknesses δ 1 and δ 2 . This method is known as the Simplified Integral Method (SIM). For a complete description of the methods, see [1]. A major point for aerodynamic solvers when computing the boundary conditions is wall roughness. Indeed, ice accretion is an iterative process during which a rough ice layer appears on the wall. The main consequence for the boundary layer is an increase in the wall friction coefficient C f and in the wall heat flux (or heat transfer coefficient h t ) magnitudes. Roughness modeling consists in taking into account its influence on the boundary layer by defining an equivalent academic well characterized roughness distribution. The similitude between the real and the academic roughness effects is based on the equivalent sand grain thickness k s [13,14].

It is a function of the geometrical and topological characteristics of the academic rough elements such as the average thickness, the typical shape (sphere, cylinder, . . . ), the distance of each other and the spatial distribution. The definition of k s and more generally the roughness influence on the boundary layer remains yet a very active research field [15,16,17,18,19,20,21,22,23,24,25]. I have not been personally involved in these studies dedicated to roughness. For a detailed review see the book of Kays & Crawford [26]. The way these methods are taken into account in the simplified integral method used in IGLOO2D (SIM2D) is described in [1].

The trajectory solver TRAJ2D is dedicated to the computation of the dispersed two-phase flow fields. For a complete description of the modeling of dispersed two-phase flows, see [27]. TRAJ2D is in charge of transporting the particles from their injection to their interaction with the wall. The deposited mass flow rate and the physical properties of the impacting particles (temperature, velocity, liquid fraction,...) will serve as inputs for the accretion solver. Two versions of TRAJ2D are proposed. Firstly, the TRAJL2D trajectory solver is based on a Lagrangian description of the problem. A detailed description of TRAJL2D is proposed in [1], itself derived from [28]. TRAJL2D allows the consideration of many physical phenomena:

• Computation of both liquid droplets and solid ice particles.

• Several drag models available which allows the computation of non spherical particles like ice particles (Sec. 4.2.1).

• Particle temperature evolution, as well as phase changes (evaporation, condensation, sublimation or melting). See Sec. 4.2.2.

• Atomization process for large droplets.

• Complex wall/particle interactions for ice particles (Sec. 4.3) and SLDs (Sec. 3.2).

• Coupling with the accretion solver (to get liquid film thickness or liquid mass fraction at the wall).

For liquid droplets and a typical LW C = 1 g.m -3 , the particle volume fraction is ≈ 10 -6 1. Therefore, no feedback from the droplets to the carrying phase (air) is performed. Secondly, the TRAJE2D trajectory solver is based on an Eulerian description of the problem. See [27] for a detailed description. The following monokinetic equations are used:

∂αp ∂t + ∇ • (α p V p ) = 0 ∂αpVp ∂t + ∇ • (α p V p ⊗ V p ) = α p V -Vp τp (5.1)
where α p is the particle volume fraction, V p the particle average velocity, V the air gas velocity and τ p the particle dynamical response time. The models proposed in TRAJE2D are less sophisticated than those in TRAJL2D. Only the dynamical description of the particles is available (no equations for temperature). For the time being, TRAJE2D can only deal with a single class for particle diameter and is restricted to App. C. Indeed, the models for complex wall/particle interactions are not available: only the full deposit mode is proposed. As the velocity in every computational cell is defined from a single averaged velocity (Eq. (5.1)), TRAJE2D is not available to deal with particle jet crossing or secondary (re-emitted) particles near a wall. To address this last point, an Eulerian multi-velocity moment model would be required [29]. That could be the subject of future works for TRAJE2D. Ice accretion is computed from the outputs provided by the aerodynamic and the dispersed twophase flow solvers. In IGLOO2D, several solvers, all included in the generic term ACCRET2D, are designed for this purpose. MESSINGER2D is a solver based on the Messinger balance [30] described in Sec. 3.1 (see in particular Fig. 3.1 and Eq. 3.1).

A few words about the general methodology used for IGLOO2D which is represented in Fig. hereafter the multistep method). From only two iterations of the loop proposed in Fig. 5.3, and except for large horn-ice shapes, accurate ice shapes can be computed (Sec. 5.3) for a moderate CPU cost with the predictor/corrector method. For a detailed description see [1]. Extensions of the classical Messinger approach were proposed. For example, Chauvin [31,32] defined a multi-layer model (MiLeS2D) composed of pure ice which may be partially melted surmounted by a runback water film. Bennani [33,34,35,36] developed a solver (ETIPS2D) dedicated to Electro-Thermal Ice Protection System (ETIPS) where the accretion solver is coupled with the heat conduction problem in the wall. Crack initiation is taken into account as well as ice shedding. Section 6 is dedicated to these extensions.

IGLOO2D assessment (App. C)

IGLOO2D has been validated through a database of more than 50 test cases [37] with different geometrical profiles, different aerodynamic configuration (velocity, temperature) and different icing condition (droplet diameter, cloud liquid water content and accretion time). This database is a mix between rime and glaze ice configurations. Most cases are covered by App. C and some configurations with SLDs have been tested as well. For rime ice shapes (Fig. 5.4), the comparison with the experimental results are very accurate [37]. Indeed, in this case, the computed ice shape derives directly from the calculation of the collection efficiency, which does not require any particular difficulty in the framework of App. C (with simple geometrical geometries). For glaze ice shapes with large horns, the results are more mitigated. To illustrate this, the two techniques described above, namely the robust but less accurate predictor/corrector technique on coarse structured grids and the accurate but more expensive multistep approach on finer unstructured grids, are tested in [1]. They are compared with the results from the NASA's icing code LEWICE. All the numerical simulations are compared to the experimental ice shapes obtained in the IRT of NASA's Glenn Research Center. The results are presented in Fig. 5.5 and the input conditions are given in Tab. 2 from [1]. The ice shapes obtained with LEWICE and the predictor/corrector technique for IGLOO2D are quite similar except near the leading edge where the ice thickness obtained by IGLOO2D is larger than the one obtained with LEWICE. This may be explained by the inaccuracy in the computation of the heat transfer coefficient near the leading edge. For glaze ice shapes with larges horns (Figs. 5.5c, 5.5e and 5.5f), even if the results are not yet very satisfactory, the multistep approach improves the horns' angular position. Regarding the horn height and width, there is still room for improvement (Figs. 5.5b and 5.5c for instance). Generally speaking, the multistep approach allows to improve the accuracy of the computed ice shapes. However, and this is the main point to date, there is no best practice to adjust the local time step ∆t/N of the multistep method as a function of geometry (grid, local curvature) or aero-icing conditions. Here, N is the total number of time steps and ∆t the total accretion time. Thus, the computed ice shapes still often depend on N .

Conclusions and perspectives

I have started the re-engineering of the icing 2D tools since 2011. The first official release was delivered to the industrial partners in 2017. It contains all the models at the state of the art for App. C, O (SLD) and D (ice crystals and mixed phase). The code architecture was improved too. However, some hot spots have not been solved yet which include:

• Automatic meshing which is still limited to simple configurations: multiple closed profiles or simplified vane rows.

• Best practices to define the number of time steps necessary to obtain an accurate ice shape independent of the grid scale effects. Although it seems to be more precise and promising, we have no hindsight on the use of the multistep approach.

• Gaps in modelling for boundary layers on rough walls and the associated heat transfer coefficient near the stagnation point. These concerns are challenging since they are of first order to estimate the final ice shape.

• The coupling with ice protection system (IPS). This point is addressed in Sec. 6 and is the next key objective for IGLOO2D.

Most existing industrialized codes propose a stationary approach where the different solvers are chained in a sequential manner: meshing, aerodynamic flow field computation, particle trajectory to finish with accretion solver and the Messinger balance. The next generation of icing codes will have to be designed with an unsteady approach necessary to simulate phenomena such as crack propagation inside an ice block or the periodic ice shedding observed with IPS in deicing mode.

Part II Models and tools dedicated to electro-thermal ice protection systems (ETIPS) were developed in the framework of two PhD studies. Firstly, Rémi Chauvin [1] proposed a three layer approach allowing to model in an unsteady way ice accretion and runback. This approach proved to be very efficient for modeling the melting of the ice layer and the derived runback ice. I supervised the PhD studies of R. Chauvin with P. Villedieu. At the same time, Lokman Bennani [2] proposed a modeling of electro-thermal ice protection systems. The mechanical behaviour of ice and fracture was studied as well as ice shedding. My involvement in his PhD works led to the writing of a common article [3]. Regarding this topic, 2 papers [3,4] and 1 conference paper [5] were written. Note that I will not talk about electro-mechanical ice protection systems (EMIPS). In the rest of the document, the term IPS will implicitly refer to electro-thermal ice protection systems.

Ice protection system modeling

Introduction

The previous chapters have shown the risks encountered by the aircrafts due to the icing phenomenon. Be it for the App. C (Sec. 5), supercooled large droplets (SLD, Sec. 3.2) or ice crystal icing (Sec. 4), aircraft manufacturers must therefore comply with certifications and regulations regarding flight safety in icing conditions. In order to achieve that goal, several ice protection technologies may be used [4]. A commonly employed system for large airliners is the bleed-air anti-icing system. However, in the context of more electrical aircraft and reduction of fuel consumption, alternative systems are currently under consideration for large airliners. The main idea is to use a system which would be able to function in de-icing mode. That is to say that a reasonable amount of ice is allowed to build up. Periodically, the system is activated and the ice is shed from the surface or melted. Such a system would have an unsteady cyclic operating mode, enabling a more optimized use of energy.

The classical tools based on the Messinger approach (Sec. 5) are not well adapted for IPS modeling for several reasons [4]:

1. The Messinger balance is a steady state model unsuitable for the deicing mode where the phenomena are intrinsically unsteady (like ice shedding for instance).

2. A major underlying assumption of the Messinger model is an uniform temperature (T w = 273.15 K) when both liquid and solid water is present on the wall. Thus, in the presence of several superposed layers of liquid and solid water, the temperature does not vary in the ice block in the direction normal to the surface. Therefore, heat fluxes are not taken into account.

For instance, the Messinger model is not able to capture the melting dynamics of an ice block when an ice protection system is activated.

3. Another major assumption of the Messinger balance is that locally, all the remaining liquid water which has not frozen or evaporated just runs back downstream. Liquid film dynamics is not taken into account.

Regarding the last point (point 3), many authors propose to account for the dynamics of the liquid film [6,7,8]. Although more sophisticated approaches have been developed [9], lubrication theory is widely used and sufficient for most applications [10,11,12]. To address points 1 and 2, Myers proposes to extend Messinger model by using a multi-layered approach (one layer of ice and two layers of liquid water) [13,14]. In each layer, a linear temperature profile is assumed in the direction normal to the surface. However, a quasi-steady hypothesis is supposed like in [15,12,16]. Later, Myers extended this idea to unsteady one dimensional cases by assuming a cubic temperature profile in each layer [17]. However, the total thickness (liquid+solid water) is fixed and evaporation or ice build up due to impacting droplets is not taken into account. Moreover, a major point is that the problem of creation/vanishing of a layer is not addressed. The purpose of the software developed by R. Chauvin in his PhD works [1] was to extend the previously mentioned research to take into account all unsteady phenomena encountered during ice build-up on a surface (eventually heated by a system) such as: water runback, dynamic phase change interfaces, unsteady temperature fields, creation/vanishing of a layer. This sofware called MiLeS2D is described in [1,4] and is the subject of Sec. 6.2. MiLeS2D is a more advanced alternative to MESSINGER2D for ACCRET2D (Fig. 6.1b).

At the same time, during his PhD works [2], L. Bennani developed a suite of tools dedicated to ice protection systems. The solver ETIPS2D is in charge of solving heat equation in the electro-thermal protection system. The latter consists of a metal thickness (the wall at the leading edge) in which heating mats are integrated to form a sandwich structure. The solvers MiLeS2D (or MESSINGER2D) and ETIPS2D are coupled with a Schwarz algorithm (Fig. 6.1a). These solvers are the subject of Sec. 6.3. Thus, the heat provided by the mats diffuses into the wall (ETIPS2D) which allows the ice to melt locally (MiLeS2D), allowing initiation and fracture propagation within the accreted ice block. The solver CRACK2D is in charge of ice fracturing and ice shedding (Sec. 6.4). Figure 6.1b outlines the core solvers of IGLOO2D, the ones dedicated to IPS and their coupling.

A triple layer based software as a substitute for the Messinger approach: MiLeS2D.

This section is strongly inspired by the PhD work of R. Chauvin [1] and from the paper [4]. Regarding MESSINGER2D, the solution of the thermodynamic problem is given by a three-step approach, where each of the regimes (rime ice, running wet and glaze ice) are tested (Sec. 5.2). As far as MiLeS2D is concerned, six distinct modes (labeled 1-6 in Fig. 6.2a) are considered [1,4]: 1. Full evaporative: the whole mass of impacting droplets is evaporated (for example, due to heat provided by an ice protection system).

2. Running wet: only a liquid water film is present. Under the action of the aerodynamic forces, the liquid film runs back along the surface.

3. Rime accretion: the droplets freeze almost instantaneously leading to ice build up with no liquid water.

4. Glaze accretion: the droplets freeze, but at a slower rate than in the rime case. Therefore, a running liquid water film is present on top of the ice layer.

5. Rime accretion with melting at the surface (due to heat provided by an ice protection system for example).

6. Glaze accretion with melting at the surface (due to heat provided by an ice protection system for example).

Physical modeling and equations

In each mode previously described, one can distinguish the presence (or abscence) of one or more layers. The term "three-layer model" comes from the fact that in the most general case one has to deal with three layers (Fig. 6.2b):

1. Running film: this is the water layer directly in contact with the airflow. It is sheared by the aerodynamic flow and runs back downstream.

2. Ice: this layer represents the ice that may build up in icing conditions.

3. Static film: when a system is activated, a part of the ice layer may melt. This has the effect of forming a static film layer (it is not in contact with the external aerodynamic flow).

Each of these three layers is described by a set of equations (see hereinbelow). However, to describe the six modes available in MiLeS2D (Fig. 6.2a), these three layers have to be combined to each other, leading to boundary conditions and source terms added to the set of equations and specific to the considered mode. Firstly, we propose to describe the equations for each layer. Then, the relevant boundary conditions and source terms for the coupling between layers are presented. Only the main ideas will be presented. For a complete description, see our dedicated paper [4] and the PhD thesis from R. Chauvin [1].

Equations for the running film layer

The running film is supposed to be very thin so that the classical lubrication theory can be applied [18]. The mass and energy equations can be written:

∂ρ w h f ∂t + ∂ρ w h f vx ∂x = Γ f (6.1a) ∂ρ w c w h f Tf ∂t + ∂ρ w c w h f vx Tf ∂x = Φf (6.1b)
where h f is the liquid film thickness, vx the mean velocity, Γ f (resp. Φf ) the mass (resp. energy) transfer term (detailed later and specific to each mode). The liquid film is assumed to be very thin and of uniform temperature on the normal direction. Therefore, in the case where the liquid film is (locally) running over ice, the previous energy equation degenerates to Tf = T m where T m is the melting temperature (T m = 273.15 K).

A TRIPLE LAYER BASED SOFTWARE AS A SUBSTITUTE FOR THE MESSINGER APPROACH: MILE

Heat conduction in the ice layer and melted film layer

Both the ice and melted film layers are assumed motionless. The mass and energy equations are respectively given by:

∂ρ k h k ∂t = Γ k (6.2)
and

∂ρ k c k T k ∂t = λ k ∂ 2 T k ∂z 2 for z ∈]a k (t), b k (t)[ (6.3a) f a T k , ∂T k ∂z = 0 for z = a k (t) (6.3b) f b T k , ∂T k ∂z = 0 for z = b k (t) (6.3c)
where the subscript k is used to denote the two layers (k = i for the ice layer and k = s for the static melted film layer). Heat transfer in the tangential direction is neglected, which explains the shape of Eq. (6.3a). Figure 6.3 shows the gradients ∂T ∂x and ∂T ∂y for a case where an internal heating is applied to a specific region of an ice block (see [4], Appendix A). The x component of the temperature gradient (Fig. 6.3a) is non zero only in a confined region at the edge of the heated area, whereas variations in the y direction (Fig. 6.3b) are much more present. This is enough to build a 1D model that retains the relevant physics. f a and f b (Eqs. (6.3b) and (6.3c)) denote in a generic way the boundary conditions at the moving boundaries of each layer, as explained in the next paragraph. This allows to combine the layers to represent each mode.

The six modes: boundary conditions, source terms and coupling between layers

The equations for each of the six modes are derived by combining those dedicated to the three layers (Eqs. (6.1), (6.2) and (6.3)) by adapting the corresponding boundary conditions (f a and f b ) and source terms (Γ f , Γ i , Γ s and Φf ). The complete description of the equations is available in [1]. A more synthetic description is also proposed in [4]. Note that in the modes with a liquid running film layer, the hypothesis of a constant temperature only is too restrictive since the computation of the temperature gradient in the running liquid film is then not possible. This is a problem to compute the rate of melting/freezing. So as to overcome this limitation, the temperature gradient is also assumed constant (and different from zero) as if the temperature profile was in fact linear [4].

Numerical method

Concerning spatial discretization, an original idea is proposed in [1,4]:

• Regarding x-wise direction, Eqs. (6.1) are discretized with a finite volume Roe scheme [19]. The method is classic and is not detailed.

• Regarding z-wise (normal) direction, the heat equation for the (eventual) ice or static film layers (Eqs. (6.3)) is discretized using a Galerkin method. This point, as well as the implicit time integration, are described in detail in [4]. A short description is proposed in this manuscript.

A Galerkine-type finite element spatial discretization in the z direction for the heat equation (Eq. (6.3))

A variational formulation of Eq. ( 6.3) is written [4]:

b k a k ρ k c k ∂T k ∂t v(z) dz - b k a k λ k ∂ 2 T k ∂z 2 v(z) dz = 0 (6.4)
where v(z) is a test function. Equation (6.4) is discretized in a subspace of finite dimension:

∀t ≥ 0 ∀z ∈] -1; 1[ T k (t, z) = n j=1 θ k,j (t)Ψ j (z) (6.5)
where (Ψ j ) n j=1 is a set of linearly independent functions and z is the dimensionless variable for z. Polynomial functions are chosen for the Ψ j [4]. Finally, the following matrix expression, whose main unknowns are the weights θ k , is obtained:

ρ k c k b k -a k 2 A θk - 1 2 ȧk + ḃk Π + ḃk -ȧk Γ θ k + 2λ k b k -a k Λθ k = f k (θ k ) (6.6)
where the definitions for the matrix A, Π, Γ and Λ and the vector f k are detailed in [1,4].

Implicit time marching algorithm

A general implicit Euler time scheme suitable for all modes is proposed in [1,4]:

M md U n+1 k -U n k ∆t + G md U n+1 L , U n+1 k , U n+1 R = S md U n+1 k (6.7)
where the matrix M md , the flux G md and the source term S md depend on the mode md. To be more precise, the flux G md derives from the spatial discretization of Eq. (6.1) (if applicable) with the Roe scheme. k is the current cell whereas L and R respectively refer to the left and right neighbors of the cell k. The source term S md is a local function (i.e. depending only on the local variables of the current cell) which models heat transfer in the normal direction with phase change. It derives from the previously described Galerkin method (Eqs. (6.2) and (6.3)). The dimension of the unknown vector U n+1 k itself depends on:

• the mode md. For instance, in the glaze mode, U n+1 k = (h f , h i , θ i,1 , . . . , θ i,ng ) where ng is the number of degrees of freedom of the Galerkin method. In the running wet mode, U n+1 k = h f , Tf .

• the spatial position k. Indeed, two different modes can be observed for two successive cells. For instance, when switching from the running wet to the full evaporative mode as the whole liquid film evaporated.

• the time n. As the three layer approach allows to take into account in an unsteady way ice accretion and runback, the mode at a give place can change between two successive times.

A fixed point algorithm is proposed in [1,4] for the implicitation of time integration of Eq. (6.7). This allows to address two major challenges. Firstly, all the cells of the grid are coupled with their left and right neighbouring cells through the running liquid film. To overcome this difficulty, a cell numbering is chosen for iterating over the cells. The idea is to run through the cells in the direction of propagation of information. It is given by the direction of shearing of the film, which determines the main direction of propagation of the film downstream. Regarding cell k, this allows to restrict the implicit formulation to the term U n+1 k itself and either to U n+1 L or U n+1 R depending on the "upwind" direction. Secondly, the implicit procedure has to manage the appearance or disappearance of a layer and the switching between modes. In a given mode md, solving Eq. (6.7) may yield a solution U n+1 k which is incompatible with this mode. To address this difficulty (which appears to be the hot spot of the method), a kind of trial process for the different modes is proposed (Fig. 5 from [4]). This very technical point is described in detail in [1,4].

Assessment and basic computations

In [4], the Stephan problem [20] was simulated to validate the Galerkine-type method. It is shown that with 8 degrees of freedom (i.e. n = 8 in Eq. (6.5)), the results computed with the Galerkine method are in good agreement with the analytical solution. It is also shown that in steady icing conditions, MiLeS2D is able to predict the same results as Messinger's method, which confirms that MiLeS2D has been really designed for unsteady configurations with a heated leading edge. As such, a test case with delayed activation of the electro-thermal protection system and runback ice build-up is proposed in [4]. After some ice had the time to build up (20 s), an internal heating is activated to eliminate it. Such a case cannot be simulated using Messinger's approach. Some results at different times are represented in Fig. 6.4a. During the unsteady melting process, small regions of melted water are released and flow downstream (Fig. 6.4a). The mechanism of the formation of small ice ridges (drawn from [4]) can be explained as follows (Fig. 6.4a):

• t = 20.03 s: in a given cell at the end of the ice block, once the melting front has propagated through the whole thickness of ice, the local mode switches to "running wet".

• t = 20.08 s: the liquid water comes into contact with the external flow field. The liquid film runs back downstream under the effect of the aerodynamic forces. • t = 20.21 s: however, as the liquid leaves the protected area, it is no longer heated enough. It starts to freeze under the effect of the cold external conditions.

• t = 20.3 s: therefore, a small ice ridge starts build up outside of the heated area.

The final formation of small ice ridges is shown in Fig. 6.4b at the boundaries of the electro-thermal protection system. This test case embodies the advantages of MiLeS2D compared to MESSINGER2D, namely the capability of modeling the unsteady activation of an electro-thermal ice protection system (notably in de-icing mode).

A model for electro-thermal ice protection system

This section is based on the PhD works of R. Chauvin [1] and L. Bennani [2]. The solver ETIPS2D dedicated to the solution of the heat equation inside the metallic wall at the leading edge (Sec. 6.3.1) is coupled with the accretion solver MESSINGER2D (Sec. 5.2) or MiLeS2D (Sec. 6.2) through a Schwarz algorithm (Sec. 6.3.2). A paper is being finalized on this subject. The generic term for electro-thermal ice protection system and its modeling refer to the coupling between ETIPS2D and MiLeS2D (or MESSINGER2D).

Thermal modeling: ETIPS2D

The electro-thermal ice protection architecture [2] is composed of several layers of different materials (Fig. 6.5b). Each layer can be defined by its own thickness and material parameters. Regarding the chordwise distribution (Fig. 6.5a), several zones associated with the location of the heating elements are defined. A structured orthogonal grid is generated (Fig. 6.5c) from which the heat equation with source terms is solved [2]:

ρc p,K |K| T n+1 K -T n K ∆t = Q src,K |K| + E∈∂K (Λ∇T ) • n f |E| (6.8)
where Q src,K stands for the source terms provided by the heater mats. The thermal conductivity matrix Λ takes into account the different thermal conductivities along the chordwise and normal directions. The way Eq. (6.8) is discretized as well as the basic validation test cases are described in [2] in detail.

Coupling between MiLeS2D and ETIPS2D with a Schwarz algorithms

Equations for the solvers ETIPS2D (Sec. 6.3.1) on one side and MESSINGER2D/MiLeS2D (Sec. 6.2) on the other are numerically solved with a time implicit integration. The two solvers (ETIPS2D and MiLeS2D for instance) have to be thermally coupled to their common interface, i.e. the wall. The coupling method is an iterative process that ensures, at convergence, both the continuity of the temperature field and the temperature gradient at the wall. To do this, a Schwarz method is used. It is based on domain decomposition. Within a generally large domain, several subdomains are identified for which subproblems independent of the other domains are defined. At the interfaces between the subdomains are defined the joining conditions that allow equivalence with the initial global problem. This has been the basis for overlapping Schwarz methods [21] for which the intersection of subdomains may be non empty. This limitation due to the necessary overlapping between the subdomains was addressed by Lions [22]. Finally, Charton [23] proposed to approximate the non-local operators introduced by Lions for the joining conditions between subdomains by local operators. This has made it possible the straightforward integration of these methods into computational codes. The strength of non-overlapping Schwarz methods in modeling electro-thermal ice protection systems is to consider the numerical solutions of ETIPS2S and MiLeS2D independently, and to take into account the communication conditions at their common interface in order to ensure, at convergence, the continuity of temperature and heat flux. The adaptation of Schwarz methods to our ETIPS configurations is described in details in [1]. It is based on the following 1D general theoretical problem:

                                         ρ 1 c 1 ∂T 1 ∂t (t, x) = ∂ ∂x λ 1 ∂T 1 ∂x (t, x) ∀(t, x) ∈ R + × Ω 1 ρ 2 c 2 ∂T 2 ∂t (t, x) = ∂ ∂x λ 2 ∂T 2 ∂x (t, x) ∀(t, x) ∈ R + × Ω 2 λ 1 ∂T 1 ∂x (t, -l 1 ) = f 1 (T 1 (t, -l 1 )) ∀t ∈ R + λ 2 ∂T 2 ∂x (t, l 2 ) = -f 2 (T 2 (t, l 2 )) ∀t ∈ R + T 1 (t, 0) = T 2 (t, 0) ∀t ∈ R + λ 1 ∂T 1 ∂x (t, 0) = λ 2 ∂T 2 ∂x (t, 0) ∀t ∈ R + (6.9a) (6.9b) (6.9c) (6.9d) (6.9e) (6.9f)
where the two geometrical domains Ω i , the unknowns T i and the physical properties ρ i , c i (assumed here to be independent of x) and λ i are defined in Fig. 6.6. Any functions f 1 and f 2 are possible. The only restriction for the problem to be properly addressed is that f 1 and f 2 are increasing functions [1].

x l2 0 -l1 T1 (t, x) T2 (t, x) ρ1(x), c1(x), λ1(x) ρ2(x), c2(x), λ2(x) Ω1 = ]-l1; 0[ Ω2 = ]0; l2[
The time discretization of Eq. (6.9) can be written:

                                       -∆ T n+1 1 + µ 2 1 T n+1 1 = µ 2 1 T n 1 ∀x ∈ Ω 1 -∆ T n+1 2 + µ 2 2 T n+1 2 = µ 2 2 T n 2 ∀x ∈ Ω 2 λ 1 ∂T n+1 1 ∂x (-l 1 ) = f 1 T n+1 1 (-l 1 ) λ 2 ∂T n+1 2 ∂x (l 2 ) = -f 2 T n+1 2 (l 2 ) T n+1 1 (0) = T n+1 2 (0) λ 1 ∂T n+1 1 ∂x (0) = λ 2 ∂T n+1 2 ∂x (0) (6.10a) (6.10b) (6.10c) (6.10d) (6.10e) (6.10f)
where

T n 1 = T 1 (t n , •), T n 2 = T 2 (t n , •), µ 2 1 = ρ 1 c 1 / (∆tλ 1
) and µ 2 2 = ρ 2 c 2 / (∆tλ 2 ). The Schwarz method for the solution of the coupled problem (6.10) consists in introducing the following mathematical series

T (k) 1 , T (k) 2 verifying:        T (k+1) 1
solution of Eqs. (6.10a) and (6.10c)

λ 1 ∂T (k+1) 1 ∂x (0) = λ 2 ∂T (k) 2 ∂x (0) + ω 1 T (k) 2 (0) -T (k+1) 1 (0) (6.11a) (6.11b)        T (k+1) 2
solution of Eqs. (6.10b) and (6.10d)

-λ 2 ∂T (k+1) 2 ∂x (0) = -λ 1 ∂T (k+1) 1 ∂x (0) + ω 2 T (k+1) 1 (0) -T (k+1) 2 (0) (6.12a) (6.12b)
where the red terms are the unknown variables to be computed at each step. It is demonstrated in [1] that the series

T (k) 1 , T (k) 2 converge to T n+1 1 , T n+1 2
for the unsteady problem with linear functions f 1 and f 2 and for the steady problem with any boundary conditions (f 1 and f 2 not necessarily linear). Optimal coefficients ω 1 and ω 2 are proposed to improve the speed of convergence of the Schwarz method (decrease of the number of iterations (k) in Eqs. (6.11) and (6.12)). For the general unsteady problem with any formulation for f 1 and f 2 , no theoretical demonstration is proposed. The convergence criteria for Eqs. (6.11) and (6.12) as well as the optimal coefficients ω 1 and ω 2 determined from the two previous simplified configurations (unsteady problem with linear functions f 1 and f 2 and steady problem with f 1 and f 2 not necessarily linear) are applied successfully in the general case [1]. Figure 6.7 shows the convergence rates obtained for the steady version of the problem (6.10). f 1 is a linear function of T which represents an interface term with the outer boundary layer through a heat transfer coefficient. f 2 deals with an evaporative term and is a non-linear function of the temperature T . In Fig. 6.7a, non-optimal ω 1 and ω 2 coefficients are applied in Eqs. (6.11) and (6.12), whereas the optimal formulation is used in Fig. 6.7b. With the optimized version for ω 1 and ω 2 , the convergence of the Schwarz algorithm is significantly increased. The extension to 2D configurations is done straightforwardly by choosing the coupling coefficients calibrated for 1D problems [1].

Mechanical modeling and ice shedding: CRACK2D

This section is based on the PhD work of L. Bennani [2] and our common paper [3]. A scenario for ice shedding when an electro-thermal ice protection system operates in de-icing mode is proposed. It is represented in Fig. 6.8a. The heater labeled C (parting strip) remains activated during the whole The whole length is melted (L f = L t ) and ice no longer adheres to the surface.

Adhesive interfacial debonding : part of the length L f = x%L t is melted and the adhesion forces that maintain ice on the surface are no longer strong enough.

Cohesive interfacial debonding : part of the length L f = x%L t is melted, ice can still adhere, but a crack may nucleate due to stress concentration and propagate along the interface.

Bulk failure : part of the length L f = x%L t is melted, ice can still adhere, but a crack may nucleate due to stress concentration and propagate inside the bulk of the ice block, therefore tearing off a part of the ice.

Ice shedding is due to an interplay of all or part of the previous mechanisms, as the experiments of Wei et al. seem to suggest [24]. Table 6.1: Ice shedding: basic mechanisms. Reprinted from [2].

cyclically so as to remove the ice which has accumulated. The activation of the heaters leads to the melting of the ice in contact with the wall and the release of a liquid water film. The role of this liquid film in ice shedding phenomenon is illustrated in Fig. 6.8b and 6.8c. The extension of the thin liquid film is denoted L f compared to the total adhesion length L t of the ice block on the wall (Fig. 6.8b).

The pressure at the contact point P c between the external flow and the film (Fig. 6.8c) is entirely redistributed by the static film over the length L f . Just after the stagnation point P c , the flow is accelerated with a pressure decrease at the same time. Thus a pressure gap is generated which favors a lifting force (Fig. 6.8c). In addition to this force, tangential viscous forces are also present, which initiates a momentum at the origin of ice shedding. The mechanism underlying ice shedding, particularly crack nucleation and propagation in the core of the ice block or at the interface between the wall and the ice block, is a combination of several more elementary mechanisms. They are described by L. Bennani in his PhD work (see Tab. 6.1 extracted from [2]). We propose to describe the brittle (bulk) failure which is the subject of a common paper [3].

The brittle (bulk) failure

This mechanism is described in detail in [2,3]. The model is constructed by adapting the approach adopted by Miehe et al. [25] to a stationary case. The idea is to introduce a parameter d that characterizes the local state of damage/fracture [26] in the ice block. The starting point is a principle of conservation of energy [27]. When a solid is deformed by action of external forces, it internally stores elastic deformation energy. If, locally, this energy exceeds a certain critical energy then it will cause an increase in crack surface [3]. Therefore, when external forces are applied, the change in energy is equal to the work produced by those forces:

δE crack + δE el = Ω f vol δu dV + Γ 1 f surf δu dΓ (6.13)
where E crack and E el are respectively the total crack and elastic energies. f vol is a volume force and f surf a surface force applied on Γ 1 ⊂ ∂Ω. The displacement field is denoted u. If the elastic strain (resp. fracture) energy per unit volume is denoted ψ( , d) (resp. φ(d, ∇d)), then: d) dV (6.14a)

E el = Ω ψ( ,
E crack = Ω φ(d, ∇d) dV (6.14b)
with = 1/2 ∇u + T (∇u) the strain. The crack energy is obtained by using a regularized crack energy functional given by Bourdin et al. [28]:

E crack = Ω φ(d, ∇d) dV = Ω g c 1 2l d 2 + l 2 ∇d • ∇d dV (6.15)
where g c is the crack energy release rate and the damage variable d lies between 0 and 1 (d(x) = 0 corresponding to an undamaged state and d(x) = 1 to a fractured state) [3]. The term 1 2l d 2 + l 2 ∇d • ∇d can be interpreted as a regularized approximation of the Dirac δ-function related to the crack surface. l is an adjustable parameter that controls the width of the regularized crack. From variational arguments (see [2] for a complete detail of the computations), we obtain:

σ = ∂ψ ∂ in Ω (6.16a) g c l d -g c l∆d = - ∂ψ ∂d in Ω (6.16b) ∇d • n = 0 on ∂Ω (6.16c)
The classical relationship ψ( ) = λ 2 tr( ) 2 + µtr( 2) valid for a linear homogeneous isotropic elastic material cannot be used as is. It is considered that in the case of brittle fracture, only the tensile energy interacts with crack propagation. Thus the elastic energy is split into purely tensile and compressive parts:

ψ( , d) = f (d)ψ + 0 ( ) + ψ - 0 ( ) (6.17)
where ψ + 0 and ψ - 0 represent respectively the undamaged tensile and compressive elastic energies. Only tensile energy interacts with crack propagation. Thus only ψ + 0 ( ), the tensile part, is multiplied by a function f (d). It represents the degradation of tensile energy due to crack formation. The choice of the function f , the way the irreversible aspect of crack propagation is preserved, the complete set of equations, the iterative algorithm to solve it as well as the finite element method used to discretized the problem, are detailed in [2,3]. The main advantage of this method is that it does not require an initial crack. However [2], the main drawback is the need of a refined mesh. Especially in the case of imposed boundary load, the crack tends to spread if the mesh is not refined enough.

Basic validation of the model for bulk failure is done in [2,3]. To illustrate the mechanism, two more realistic test cases are proposed: a flate plate and a NACA0012 airfoil in an airflow and on the surface of which ice initially accreted. The experiment consists in varying L f (the extension of the thin liquid film, see Fig. 6.8b) until fracture occurs, starting with L f = 0. The configuration is proposed in Fig. 6.9a for the ice attached to the flate plate. In Fig. 6.9b, the aerodynamic forces are strong enough as well as the ratio L f /L t to initiate and propagate a crack. On the other hand, in Fig. 6.9c, the ratio L f /L t is reduced (keeping the same magnitude for the aerodynamic forces) so that fracture does not occur. The same conclusions are drawn for the ice attached on the NACA0012 profile [3]. Figure 6.9: Fracturing. Flat plate test case. Reprinted from [3].

The brittle (bulk) failure is a possible detachment mechanism, modeled here by using damage mechanics. Leaving aside the classic limitations of the proposed model, namely for instance that the ice shape may grow in time or that the surface tension forces between the thin liquid film and the ice core may play a role, there are assumptions that have been overlooked but are of the very highest order:

• The mechanical behavior of atmospheric ice derives from empirical relations. Very few studies on the subject exist. Most studies are interested in the tensile or compressive strength but do not provide many information on mechanical characteristics in the form of well defined laws.

• Emphasis has been put on bulk failure without taking into account adhesion forces. Most likely, ice shedding is a combination of these two phenomena (Tab. 6.1).

Both points (mechanical properties of atmospheric ice and adhesive debonding) were not a priority in my research activities. I suggest the interested reader to consider the work of L. Bennani [2,29].

Conclusions and perspectives

Modeling of electro-thermal ice protection system at ONERA has made a great leap forward thanks to the PhD works of R. Chauvin and L. Bennani. The outcome of their collaboration consists in a deicing application illustrated in Fig. 6.10 where the solvers IGLOO2D, MiLeS2D, ETIPS2D and CRACK2D are used together. As detailed previously, ETIPS2D and MiLeS2D are coupled through a Schwarz algorithm. A 80 s cycle is observed for the ignition and extension of the heating mats (Fig. 6.10a). The heating mat D (the so called parting strip) is never switched off. The next heating mats are cyclically activated according to the time sequence detailed in Fig. 6.10a. For the first 80 s, the air is dry. The droplets are injected at t = 80 s. Ice shapes are computed with the multistep method (Sec. 5.2) and the outer quantities are computed every 10 s. Ice shedding is estimated every 1 s with CRACK2D. From t = 80 s to t = 100 s, i.e. 20 s after the beginning of ice accretion (Fig. 6.10b), only the sparting string is activated. This ensures that the leading edge is warm enough to prevent ice from accumulating. A runback liquid film appears and freezes on an unprotected area (heating mat E, Fig. 6.10b). Ten seconds later (at t = 110.5 s, Fig. 6.10c), the thickness of the ice block continued to increase (from t = 100 s to 120 s). However, as the heating mat E has been activated for 10 s, a thin film of liquid water is formed between the wall and the ice block. Then, the ice block is removed by aerodynamic forces (not shown in Fig. 6.10). After 120 s, the heating mats C and E are switched off and an ice block begins to form above the heating mat E. After 160 s (equivalent to 2 complete cycles), the configuration of Fig. 6.10d is obtained. The cycles follow one another and the ice blocks are released as the heating mats are switched on. However, it must be made sure that at the furthest downstream heating mat, the residual liquid water film after ice shedding is completely evaporated so as not to create refreezing ice shapes downstream. Indeed, these shapes are very damaging to the aerodynamic performance of an aircraft. The configuration shown in Fig. 6.10 is based on an innovative technique which is the first to be able to describe the operating mechanism of a deicing system in an unsteady formulation. Ice sheddig prediction is not based on an empirical correlation, but on the computation of crack initiation and propagation (CRACK2D). Future research proposals, some of which have already been initiated, may be:

• The mechanical properties of accreted ice are still not well known. Dedicated studies have been initiated in the STORM (regarding adhesion), CORICE and upcoming TRICEPS projects.

• The only possible topology with MiLeS2D is limited to the superposition of three ordered layers (static film, ice block and runback liquid film). This is well suited for deicing configurations with wall heating where the melting front is located at the wall. However, for the time being, it is not possible to take into account ice porosity with a movement of the liquid front from the wall towards outside (imbibition phenomenon experienced for ice crystal icing). This may be further developed in the MUSIC-haic project.

• We are beginning to focus on the modeling of electro-mechanical ice protection systems (EMIPS). Postdoc internships as well as a PhD thesis project dedicated to this subject are in progress.

• ONERA is beginning to get involved in passive protection systems. Combined with active systems, they can improve the system efficiency.

Introduction

The motion and stability of liquid thin films and droplets which wet a solid substrate are present in a lot of natural and industrial processes and have been the object of a lot of research studies for several decades [6,7,8]. In the context of deicing, when a thermal protection system is activated, the 101 supercooled water droplets impacting an aircraft surface do not freeze instantaneously and can coalesce and form a thin liquid film as a result of aerodynamic forces. Experimental studies show that this liquid film is not always stable and can split into rivulets that may refreeze on unprotected surfaces [9]. The modeling of rivulet flows and the accurate prediction of wet and dry surfaces is important since it has a direct influence on the wall heat and mass fluxes such as evaporation or exchanges with the boundary layer.

To do this, different approaches are proposed. A macroscopic approach, named Minimum Total Energy (MTE) criteria, and proposed by Al-Khalil [10] and da Silva [11], provides a wetness factor without explicitly simulating rivulets formation. This approach is well suited for industrial computations but it suffers from a lack of sound theoretical background and has not been so far much validated. The approach we use consists in averaging the Navier Stokes equations over the fluid thickness. The hypothesis of thin films with an assumed small ratio between the film thickness h and any in-plane length-scale is done (long-wave approximation). This can lead to a single non-linear fourth order partial differential equation for h also called lubrication equation, or a third order system of two equations for h and hu also called shallow water equations. One of the objectives is to reduce the order of the equation system (up to second order) so that it can be used in an 3D finite volume solver based on unstructured grids. To do this, an extended shallow water model is proposed (Sec. 7.2). The proposed model accounts for viscosity, capillary, and gravitational effects. To take into account the molecular forces in the vicinity of the contact line with possible film breaking into rivulets and partial wetting, an additional pressure term introduced by Frumkin and Derjaguin [12,13,14,15], called disjoining pressure, is added to the equation system (Sec. 7.3).

Modeling of thin liquid films with no contact line

In this section, it is assumed that there is no contact line. The models as well as their properties are presented.

Lubrication equation

The following assumptions are done:

• Incompressible fluid ρ = cste • Small to moderate Reynolds number Re = ρu 0 h 0 µ ∼ 1 or less • Negligible wall curvature effects
• Isothermal flow (constant surface tension coefficient, density, viscosity, . . . )

A parabolic velocity profile in the direction normal to the wall is assumed in the liquid film (Poiseuille flow):

u (z) ≡ τ a µ z - 1 2µ [∇p a + γ∇κ + ρg n ∇h -ρg t ] z (2h -z) (7.1)
where only the z-dependency is mentioned (Fig. 7.1). τ a and P a are respectively the shear stress of the outer air and the air pressure. γ is the surface tension coefficient between air and water, κ is the total curvature of the film free surface given by:

κ = -∇ •    ∇h 1 + ||∇h|| 2 1/2    (7.
2)

The gravity components g n and g t are given by g n = -g z and g t = g + g z e z (Fig. 7.1). Lubrication equation is then obtained:

∂h ∂t + ∇ • hu = 0 (7.3a) u = h 2 3µ 3τ a 2h + ρg t -∇ [P a + γκ + ρg n h] (7.3b)
which is a 4 th order partial differential equation. u is obtained from the normal (z component) integration of Eq. (7.1).

Shallow water like system (non conservative formulation)

The starting point is the following shallow water like system where both h and hu are considered to be independant variables:

∂h ∂t + ∇ • (hu) = 0 (7.4a) ∂ (ρhu) ∂t + ∇ • (ρhu ⊗ u) = ρhg t -h∇ [P a + γκ + ρg n h] + 3 2 τ a - 3µu h (7.4b)
where u stands for the averaged (z component integration) velocity field u. Equation (7.4b) (using Eq. (7.4a)) may be rewritten in an equivalent way as:

ρ du dt = ρg t -∇ [P a + γκ + ρg n h] + 3τ a 2h - 3µu h 2 (7.5)
or again, using Eq. (7.3b):

du dt = u -u τ ν (7.6)
where τ ν = ρh 2 3µ . For a macroscopic observation time scale T = L/u 0 , we can write:

τ ν T = h 3L Re L 1
according to the aforementioned hypothesis (7.7) and then: u ≈ u (7.8) This proves the formal equivalence between the lubrication equation (Eq. ( 7.3)) and the shallow water like system (Eq. (7.4)).

A remark about the convective term in Eq. ( 7.4b) where the assumption u ⊗ u ≈ u ⊗ u is done. Given the Poiseuille hypothesis for the velocity profile (Eq. (7.1)) necessary to close the source term 3τa 2h -3µu h 2 in Eq. ( 7.4b), it should rather be: u ⊗ u ≈ 6 5 u ⊗ u (+ other terms). However, while the velocities remain moderate (Re ∼ 1 or less), this has no impact on the phenomena we are concerned with. Our objective here is not to accurately model inertial phenomena such as rolling waves for instance.

Modeling of thin liquid films with a contact line

At this stage, the system (7.4) is not able yet to take into account film dewetting. Indeed, we have so far neglected any interaction between the interfaces when the film thickness h is close to 0. These interactions are at the origin of wetting/dewetting phenomena.

Disjoining pressure and energy

The most popular method in the literature to address the effect of these interactions is to define a new boundary condition at the liquid/gas interface for the pressure field being written:

P (x, y, z = h) = P a (x, y) + γκ(x, y) + ρg n h(x, y) -π d (h) (7.9)
which corresponds to the pressure definition given in Eq. (7.4b) with the addition of a disjoining pressure term π d (h) driven by the interactions between the solid/liquid and liquid/gas interfaces in the vicinity of the contact line. Following the Derjaguin notation [13], π d is defined as:

π d (h) = - de d dh (h) (7.10)
where e d (h) is the disjoining energy per unit surface area. Thus, the energy of the film per unit surface area is given by:

e f ilm =
where γ sl is the surface tension coefficient between the wall and water. To ensure that e d (h) is able to address partial wetting situations with both wet and dry areas, it is necessary for the film energy formulation (Eq. (7.11)) to degenerate into the one of a dry substrate when h = 0 and into the one of a wet substrate when h is large compared to the molecular action radius R (Fig. 7.2). Thus:

• For dry conditions, the only energy involved is that of a dry substrate γ sg . Therefore, from Eq. (7.11):

γ sg = e f ilm h = 0 ||∇h|| = 0 = γ sl + γ + e d (0) (7.12)
and e d (0) = γ sg -γ sl -γ = S = γ (cos θ s -1) where S is the spreading parameter (here negative) and θ s the static contact angle. In this figure, for the sake of clarity and without loss of generality, the film is assumed to be in static equilibrium on a horizontal flat plate (e k = 0, Eq. (7.11)). Reprinted from [1].

• For wet conditions such as h R, the disjoining energy e d must be 0.

For S < 0 and for partially wetting configurations, the disjoining energy e d must be a concave function [16]. We introduce the following simple expression that meets the previous boundary conditions for e d : e d (h) = S exp (-h/h ) = γ (cos θ s -1) exp (-h/h ) (

The thickness h represents the molecular forces range. It has to be calibrated for numerical purposes. Finally, the disjoining pressure is defined as:

π d (h) = - de d dh (h) = S h exp (-h/h ) (7.14)

An augmented formulation of the system

Equation (7.4b) is a 3 rd order partial differential equation. To reduce the order of the system and based on an idea from [17], a new independant unknown p = ∇h is added. An additional equation is obtained by applying the gradient operator to Eq. (7.4a). Finally, the augmented system may be written:

∂h ∂t + ∇ • (hu) = 0 (7.15a) ∂p ∂t + ∇ (p • u) = -∇ (h∇ • u) (7.15b) ∂ (ρhu) ∂t + ∇ • (ρhu ⊗ u) = ρhg t -h∇ [P a + γκ + ρg n h -π d (h)] + 3 2 τ a 1 - b 3 (h + b) - 3µu h + b (7.15c)
which is of 2 nd order and where the curvature κ is now obtained from:

κ = -∇ •    p 1 + ||p|| 2 1/2    (7.16)
with b the slip length preventing viscous dissipation divergence for vanishing film thicknesses [18].

Shallow water like system (conservative formulation)

In order to use a finite volume approach, a conservative formulation of Eq. (7.15) is proposed. The following variables ϕ h , ϕ p and ϕ hu are defined as:

ϕ h ≡ ∇e g+d (h) ≡ ∇ (e g + e d ) (h) = ρg n h -π d (h) (7.17a
)

ϕ p ≡ ∇e c (p) = γ p 1 + ||p|| 2 (7.17b) ϕ ρhu ≡ ∇e k (ρhu) = u (7.17c)
The Legendre transformations are defined as:

L h (h) ≡ hϕ h -[e g (h) + e d (h)] = 1 2 ρg n h 2 -hπ d (h) -e d (h) (7.18a)
L p (p) ≡ p ⊗ ϕ p -e c (p)I d = γ p ⊗ p 1 + ||p|| 2 -γ 1 + ||p|| 2 + γ sl I d (7.18b) L ρhu (ρhu) ≡ ρhu ⊗ ϕ ρhu -e k (ρhu)I d = ρhu ⊗ u - 1 2 ρh||u|| 2 I d (7.18c)
where I d the identity matrix. It can then be obtained:

h∇ (ρg n h -π d (h)) = ∇L h (7.19a) h∇ (γκ) = ∇ • (L p -h (∇ • ϕ p ) I d ) (7.19b)
The demonstration of Eq. (7.19a) is obvious. Regarding Eq. (7.19b), it is done in [5]. Therefore, the system (7.15) can be written in the following conservative form:

∂h ∂t + ∇ • (hu) = 0 (7.20a) ∂p ∂t + ∇ (p • u) = -∇ (h∇ • u) (7.20b) ∂ (ρhu) ∂t + ∇ • (ρhu ⊗ u) + ∇ • (L h I d + L p ) = ∇ [h (∇ • ϕ p )] + ρhg t -h∇P a + 3 2 τ a 1 - b 3 (h + b) Sext - 3µu h + b (7.20c)
where S ext stands for the source terms. A classical finite volume approach is chosen to discretize Eq. (7.20) and the numerical schemes (for both space and time) are described in [1,3]. It is shown that Eq. (7.20) may have complex eigenvalues due to the concave shape of the disjoining energy (Eq. (7.13)). Therefore, classical schemes used for hyperbolic equations cannot be implemented [3].

Associated energy equation

Multiplying the previous system (7.20) by the entropic variables (Eq. (7.17)) leads to:

< ϕ h ; ∂h ∂t > = ∂e g+d ∂t (7.21a) < ϕ h ; ∇ • (hu) > = ∇ • (e g+d u) + L h (∇ • u) (7.21b) < ϕ p ; ∂p ∂t > = ∂e c ∂t (7.21c) < ϕ p ; ∇ (p • u) > = ∇ • (e c u) -e c (∇ • u) + ∇ • [(p ⊗ ϕ p ) • u] -< u; ∇ • (p ⊗ ϕ p ) > (7.21d) < ϕ p ; ∇ (h∇ • u) > = ∇ • [h (∇ • u) ϕ p ] -h (∇ • u) (∇ • ϕ p ) (7.21e) < ϕ ρhu ; ∂ρhu ∂t > = ∂e k ∂t + 1 2 ρ||u|| 2 ∂h ∂t (7.21f) < ϕ ρhu ; ∇ • (ρhu ⊗ u) > = ∇ • (e k u) + 1 2 ρ||u|| 2 ∇ • (hu) (7.21g) < ϕ ρhu ; ∇ • (L h Id + L p ) > = ∇ • (L h u) -L h (∇ • u) + < u; ∇ • L p > (7.21h) -< ϕ ρhu ; ∇ [h (∇ • ϕ p )] > = -∇ • [h (∇ • ϕ p ) u] + h (∇ • ϕ p ) (∇ • u) (7.21i) -< ϕ ρhu ; 3µu h + b > = - 3µ h + b ||u|| 2 (7.21j) < ϕ ρhu ; S ext > = S ext • u (7.21k)
where < • ; • > denotes the canonical scalar product. By summing each of the equations of the system (7.21), we obtain the following equation for the film energy e f ilm .

∂e f ilm ∂t +∇•(e f ilm u)+∇•[(L p + L h Id) • u]+∇•[h (∇ • u) ϕ p ]-∇•[h (∇ • ϕ p ) u] = - 3µ h + b ||u|| 2 +S ext •u (7.
22) Provided that the source terms S ext of Eq. (7.20) are not driving forces, Eq. (7.22) ensures that the energy e f ilm associated to the system (7.20) is a decreasing function over time.

Model calibration

The behaviour of the film described by Eq. (7.20) is controlled by two fundamental parameters:

• h (Eq. (7.13)) which represents the molecular forces range in the vicinity of the contact line.

• b which is the slip length preventing viscous stress divergence for vanishing film thicknesses.

We proposed a 1D calibration of these two parameters h and b in [1,3] with respect to the local spatial step ∆x.

Influence of h (Eq. (7.13))

The static steady state of the spreading of a liquid film on a horizontal flate plate is studied [1,3]. Configuration is shown in Fig. 7.3a. Film thickness, static contact angle and spreading at equilibrium are denoted respectively h puddle , θ s and W s . The spatial step ∆x being set, the influence of h on the model capability to accurately predict θ s and W s is studied. If a sensitivity analysis with respect to h is carried out, gradually increasing its value, different regimes can be defined: • If h ∆x: spatial resolution is too coarse for the length scale h . The molecular forces in the vicinity of the stagnation line are not properly taken into account. Neither θ s nor W s are accurately computed (Fig. 7.3c).

• If ∆x h : Spatial resolution is optimal for the molecular forces derived from e d . Both θ s and W s are accurately computed (Fig. 7.3c).

• If ∆x h h puddle (typically h = h puddle /10): the length scale on which the disjoning energy is regularized is large compared to ∆x. Thus, the forces derived from the disjoining energy are spread out around the contact line. Spreading W s is computed accurately whereas the contact angle θ s is not (Fig. 7.3c).

• If h h puddle : the length scale on which the disjoning energy is regularized is too large and neither θ s nor W s are accurately computed (Fig. 7.3c).

As the disjoining energy e d stands for the interacting energies at the molecular scale, h has to be chosen as small as possible. Therefore, in view of what is shown previously, the condition ∆x h must be respected for an accurate estimation of both θ s and W s . However, if ∆x cannot be chosen too small due to computational time reasons, we show in [1] that the less restrictive condition h ∼ h puddle /10 with larger values for ∆x (typically in the order of h puddle ) is enough for an accurate estimation of the spreading W s only. Be careful, however, with excessive h /h puddle ratios, which can lead to unstable liquid film regimes with dewetting phenomena that have no physical justification [1,3].

Influence of b (Eq. (7.20c))

Following the same approach as Troian [19], Bertozzi [20], or Kondic [21], a liquid film injected with a constant flow rate uniformly on an inclined flat plate is considered [3]. After the fluid is released, and after initial transients, the flow develops a traveling wave solution, caracterised by the presence of a bump of thickness h bulk near the contact line, and a Nusselt flow of thickness h N u far from the contact line (Fig. 7.4a).

The objective is to study the influence of the slip length b (Eq. (7.20c)) on h N u and h bulk . The spatial step ∆x and the thickness h are set so that h ∼ h N u /10 and ∆x h . At first glance, be it for totally or partially wetting films, b has a great influence on the thicknesses h N u and h bulk (Figs. 7.4b and 7.4c). Actually, the conclusions are different whether a fully or partially wetting film is considered.

Regarding fully wetting films (θ s = 0 • ), Fig 7.5a shows that while an asymptotic value for h N u is quickly obtained for the small values of b (b ≤ h N u /10), the thickness h bulk continues to increase for even smaller values of b. This is confirmed by the relative error r b between h bulk and h N u (Fig. 7.5b). The kind of plateau for h bulk observed for b/∆x < 1 (Fig. 7.5a) does not suggest a converged configuration but rather a poor spatial discretization of the viscous term.

Regarding partially wetting films (θ s > 0 • ), both thicknesses h bulk and h N u increase when b decreases (Fig. 7.6a). An asymptotic value for both h bulk and h N u is reached when b ≤ h N u /10. This is confirmed by Fig. 7.6b where the relative error r b between h bulk and h N u is almost independent of b. Beware of the y-scales which are different between Figs. 7.5b and 7.6b. Therefore, all we need is to choose b small enough to ensure a significant stabilising viscous term in the vicinity of the contact line. Typically, b is chosen equal to 10 -9 m, which is much lower than ∆x as far as our applications are concerned.

3D model validation

Two test cases are studied: the pinching of a falling film and the transition into rivulets. Regarding numerical parameters, the conditions h /h bulk ≈ 1/10 or h /h N u ≈ 1/10 are respected. However, Table 7.1: Evolution of the spacing w between the rivulets and their width d according to the Reynolds number Re of the liquid injection and the inclination β of the flate plate. Experimental data obtained by Johnson [23] with fluid B. Reprinted from [1] the condition ∆x h is more challenging due to computational time restrictions. Indeed, the time implicit formulation of the 3D solver is not yet available. Nevertheless, even for spatially under-resolved configurations, the matching between numerical and experimental results is proven to be good. Note that for the 3D computations, the solution of Eq. (7.20b) does not ensure that ∇ ∧ p = 0, or, in other words, that p can be written as a gradient. A trick is to enforce the assumption p = ∇h. Ensure that ∇ ∧ p = 0 is one of the future improvements to be made in the code. In particular, some ideas based on the Helmholtz decomposition that any sufficiently smooth, rapidly decaying vector field in three dimensions can be resolved into the sum of an irrotational (curl-free) vector field and a solenoidal (divergence-free) vector field are in progress.

Pinching of a falling film

A liquid film is injected with a constant volume flow rate Q through a slit. The numerical simulations are compared to the experimental results obtained by B. Thoraval during his PhD work [22]. Six configurations are tested (3 different fluids with two volume flow rates) and are detailed in [1]. The plane is vertical with a falling liquid film. The following numerical parameters are chosen: b = 10 -9 m, h /h bulk ≈ 1/18 and ∆ x = ∆ y = 3h . The results with liquid water and two different volume flow rates are presented in Fig. 7.7. The input static contact angle θ s is measured experimentally and estimated between 60 • and 80 • . The pinching lengths, symbolized by the white solid line, is compared between the numerical and experimental results. For each of the two volume flow rates Q, the length measured experimentally is between the two ones obtained from the computations. The latter are carried out with the two extreme values of the confidence range for the experimental measurement of the static contact angle θ s (68 • ≤ θ s ≤ 80 • ). The same kind of results are obtained with the other two fluids [1].

Film transition into rivulets

The capability of the model to take into account film transition into rivulets is demonstrated in [1,3] from the experimental studies of Johnson [23]. A liquid film is injected uniformly onto a flate plate at different angles β. Different Reynolds numbers Re are observed. The spacing w between the rivulets and their width d are measured. The different configurations are summarized in [23]. One of them at Re = 0.52 is illustrated in Fig. 7.8. The following parameters are chosen for the computations: b = 10 -8 m, h = h N u /10 and ∆ x = ∆ y h N u . The spacing w between the rivulets and their width d are compared between the numerical and experimental results in Tab. 7.1. For the particular Reynolds number (Re = 0.52) and regardless the poor spatial resolution, a good agreement is found between computations and experiments. Similar conclusions are driven for the other Reynolds numbers. A more detailed study can be found in [1].

Conclusions and perspectives

A shallow water type model to simulate dynamics of thin liquid films with a contact line has been proposed. To do this, an additional disjoining pressure term has been introduced in the system. An augmented 2 nd order system has been proposed from which an associated energy equation has been derived. The system, written in a conservative way, is well adapted for a finite volume numerical discretization. Calibration of two fundamental parameters h and b has been carried out and finally 3D test cases have been done. The next steps are first of all of a numerical nature: 3D time implicit formulation and enforcement of irrotationality for p in 3D. Regarding modeling, hysteresis phenomena have to be investigated. The case of a sheared air flow above the film shoud be more thoroughly investigated too. A coupling with the thermal equation is needed to take into account a possible freezing of the liquid film.

The PhD works of J. Lallement [1] have led to improvements in the treatment of phenomena near the contact line. Progress has been made, especially in the mathematical formulation. However, the very mechanism of rivulets formation remains to be clarified. Simplified approaches such as linear stability analysis have made it possible to highlight some destabilizing parameters like the ratio h bulk /h N u . However, transition to rivulets seems to be a highly non-linear phenomenon. This third part of the manuscript is dedicated to my research activities which are not connected to icing. Chronologically, they are related to the first part of my career at ONERA. This part is composed of three chapters:

• Chapter 8: Numerical methods applied to the direct numerical simulation (DNS) of interfacial multiphase flows. This is a brief summary of my PhD works conducted between 2006 and 2009 at ONERA.

• Chapter 9: isogeometric analysis, which was one of my research topics during my first years at the DSNA (ex DAAA) department.

• Chapter 10: focus on the adjoint method for optimal aerodynamic design and error control.

tension coefficient. A parametric study based on the Weber number (W e ) was performed. Finally, an a priori study was carried out to derive different subgrid scale terms from the filtering of the DNS database (Sec. 8.4).

Front capturing methods and dedicated numerical methods

The numerical simulations were performed with DYJEAT [13], a two-phase incompresible Navier-Stokes solver. Incompressibility is ensured by a projection method [14,15] on a MAC Cartesian staggered grid [16]. Originally, the interface was advected by a level-set front capturing method [17,18]. Later, a Coupled Level Set and Volume Of Fluid method (CLSVOF) method [19,20] was proposed in DYJEAT. Jump conditions at the interface are addressed by a ghost fluid method [21].

The main outcomes presented in this section are reported in [4,2,1]. The context being primary atomization, the numerical methods have to account for a largely deformed interface with strong local curvatures. Thus, efforts were focused on the numerical methods dedicated to interface advection and on the processing of jump conditions at the interface (ghost fluid method [21]). Regarding interface advection, front-capturing methods (level-set [17,18], VOF [22]) and front-tracking methods (marker methods from Shin & Juric [23]) were presented. It was the opportunity for us to develop an hybrid method combining a level-set approach coupled with particles describing the interface [4,1]. The test cases used to validate the different numerical methods met different requirements to validate their capabilities in steep configurations encountered in atomization processes. These included turbulent flows, strong interfacial deformations (stretching, ligament formation, break-up), strong local curvatures, large density and viscosity ratios at the interface with high surface tension. Thus, in order to assess the capability of the different methods to adequately transport thin structures and sharp corners, the solid body rotation of the Zalesak's [24] notched circle was simulated (Fig. 8.1a). The Bell et al. [27] deformation of a circle by a single vortex was considered to assess the ability of the numerical methods to resolve thin filaments (Figs. 8.1b and 8.1c). The numerical methods as well as acronyms are detailed in [1]. The hybrid level-set method coupled with markers (called LSET-PART in Fig. 8.1) is the most accurate in terms of numerical diffusion and dissipation (Fig. 8.1). However, the latter may remain disadvantageous with high computing times. A detailed comparison between the different Eulerian and Lagrangian front capturing methods as well as the technique used for spatial discretization is proposed in [1,4].

DNS and interaction between turbulence and interfacial scales

The work presented in this section is a contribution to the understanding of the coupling between turbulence and capillary effects. The main outcomes are reported in [3,1]. A free HIT (homogeneous isotropic turbulence) was chosen to conduct this study. First, it was established that with a spatial resolution of 512 3 and a 5 th order WENO scheme for the spatial discretization of the non-linear terms, a spectral precision of our DNS is achieved for the single-phase HIT flow. The evolutions of the characteristic turbulence quantities (turbulence kinetic energy, viscous dissipation rate, Taylor and Kolmogorov scales, . . . ) are in line with the theory and with the ones obtained with a spectral code. Then, an initially flat sheet was added to the flow, and a parametric study based on the surface tension coefficient (W e ) was conducted. Density and viscosity ratios were set to 1. From a qualitative point of view, two kinds of interactions between turbulence and interface scales were observed. For large Weber numbers, the interface is torn by surface tension forces and fragmented into a large number of (a) Zalesak test case [24]. 64 2 grid. Black: initial (t = 0s) ; red: level-set [17,18] method with a WENO scheme [25] for the spatial discretization of the advection equation for the levelset Φ ; blue: VOF method [22] with PLIC reconstruction [26] of the interface ; green: level-set method with a spectral method for the spatial discretization of the advection equation for the level-set Φ [1] ; orange: hybrid level-set method coupled with particles describing the interface (LSET-PART) [4].

(b) Bell et al. test case [27]. 64 2 grid. Black: initial (t = 0s) ; red: front-tracking (with markers) method of Shin & Juric [23] ; blue: hybrid level-set method coupled with particles describing the interface (LSET-PART) [4]. [24] and Bell et al. [27]. The numerical methods used are detailed in [1]. Reprinted from [1].

droplets. The turbulent kinetic energy is then used for the droplet coalescence. A region of anisotropy and vorticity production develops in the vicinity of the interface. An energy transfer happens from the larger turbulent scales to the smaller ones. Its characteristic scales are related to the size of the droplets. When the Weber number is reduced, the droplets are much larger and do not interact with each other. At sufficiently long times, the interface then transfers energy back to the fluid. For even smaller Weber numbers, the interface remains almost flat with surface waves.

The work carried out in this section has been the first to propose DNS with a fully turbulent flow in which a largely deformed interface develops. Our work on understanding interface/turbulence interactions has been a first contribution. Indeed, the main principles of interface/turbulence interaction were identified:

• Three Weber ranges (strong Weber numbers with strong interface deformations, weak Weber numbers with surface waves developing on the interface and transition Weber numbers, Fig. 8.2a).

• Anisotropy and vorticity production zone around the interface (Fig. 8.2b).

• Energy transfer from large to small scales (Fig. 8.2c).

The study was limited to density and viscosity ratios equal to 1. Future work should include higher density and viscosity jumps.

DNS filtering and a priori analysis

The a priori analysis of the DNS database presented in Sec. 8.3 was performed. The main outcomes are reported in [5,6,1]. For the mass and momentum equations, a phase conditioned filtering was used. This was done to avoid filtering through the interface which represents a discontinuity for the stress tensor. For the level-set advection equation, a classical volume filtering was used, since the level-set function is perfectly continuous through the interface. Far from the interface, the classical results of single-phase flows are found and focus was given to the vicinity of the interface. The subgrid terms with the highest magnitude are the ones derived from the non-linear term (as well as for single-phase flows), and the ones derived from the coupling between the two phases. These last subgrid terms cannot be neglected. On the contrary, they are even predominant when the interface shows large deformations. Indeed, this is the subgrid contribution of the normal vector to the interface that controls the amplitude of the subgrid terms related to inter-phase coupling. The more the interface is disturbed, the higher the subgrid contribution of the normal vector. A scale similarity Bardina type model [28,29] was used for the subgrid terms derived from the non-linear terms of both the momentum equation (Fig. 8.3a) and the advection equation of the level-set function (Fig. 8.3b). The model is accurate, except in the vicinity of the interface where they do not embed anisotropy that develops close to the interface.

Conclusions and perspectives

Perspectives were proposed in the conclusions of my PhD thesis. Since I was involved in other research topics after my hiring at ONERA (firstly optimization, then icing activities), I have unfortunately not been able to take part in them personally.

Regarding front capturing methods, the tendency is now to use hybrid level-set/VOF or levelset/VOF/particle methods for interface transport which combine the advantages of each of the major approaches. At ONERA, D. Zuzio proposed a parallel adaptive algorithm to solve incompressible twophase flows with surface tension [30]. Moreover, a consistent mass-momentum advection algorithm is proposed in [31] where high density ratios can be taken into account consistently.

Regarding DNS of multiphase flows, the growing capabilities of computing resources have made it possible to perform more realistic numerical simulations with larger Reynolds numbers. As far as the interface/turbulence interaction is concerned, the density and viscosity ratios were increased. So as to keep a sufficient energy level in the computational domain, a forced HIT flow was used instead of a free HIT flow. This work is still being carried out by J.L. Estivalezes at ONERA.

Regarding the a priori analysis of DNS and the modeling of the subgrid scale terms, S. Vincent, with whom I worked during my PhD works, has performed similar works for a phase inversion configuration [32]. Different existing models such as Smagorinsky, Wall-Adapting Local Eddy-viscosity (WALE) model, Bardina, Mixed and Approximate Deconvolution Model (ADM) were used to account for twophase subgrid effects. These models were compared to filtered DNS results.

Modeling the subgrid scale terms of the coupling between turbulent and interfacial scales in the vicinity of the interface remains a challenging issue. Authors, such as M. Herrmann [33], performed a LES simulation while keeping a DNS space resolution close to the interface. The DNS field at the interface is then filtered. The reconstruction of the DNS velocity field in the vicinity of the interface remains a challenge. Over the last few years, new challenges have been emerging in the CFD community. With the rise of high-performance computing (HPC) and the reorganization of the simulation tools into effective CFD platforms whose different modules form bricks coupled to each other, ONERA has also had to adapt and develop modular high performance codes. New challenges also include meshing strategy for complex geometries whose accuracy must be comparable with that of high order numerical solvers. It is in this context that I developed a CFD code based on isogeometric analysis. It was one of my research topics during my first years at the DSNA (ex DAAA) department at ONERA. The code was written from scratch and limited to perfect fluid flows (Euler equations). A description of the isogeometric method as well as a basic assessment is proposed in [1] (conference paper).

Introduction

Isogeometric analysis is born from the need to adapt Computer-Aided Design (CAD) to Computational Fluid Dynamics (CFD) and computational mechanics. Indeed, CAD files must be translated into meshes which are analysis-suitable geometries for CFD (finite elements for instance). In the context of aeronautical industry, up to 80 % of overall analysis time is devoted to mesh generation [2]. Communication between the analysis-suitable mesh and CAD can be so complicated that some geometric operations like adaptive refinement are still dedicated to academic applications [3]. Therefore, the simultaneous treatment of CAD and CFD appears to be the solution with the use of common basis functions on which both geometry and aerodynamics fields are projected.

Another difficulty results in the approximated mesh built from the exact CAD description. Resulting errors can spread over CFD results. Indeed, Bassi [4] shows that a geometrical representation which takes into account the curvature of the boundary is mandatory in order to obtain meaningful numerical solutions. In the context of discontinuous finite elements with high order schemes, accurate numerical 131 results are challenging with geometrically linear elements even on highly refined grids [4]. Therefore, the classical piecewise linear approximation of curved boundaries is not suitable in the context of high order methods. The main solution is to focus on the CAD geometric model to be used directly as an analysis model. This concept is referred to as isogeometric analysis and was first introduced by Hughes [2]. Isogeometric analysis is an extension of classical finite elements analysis [5]. For a complete review of the practical use of the method, see [1].

My objective was to illustrate the main advantages of the isogeometric analysis for typical aeronautical CFD problems. Indeed, complex spatial configurations with curved boundary conditions are encountered both in external (wings with flaps and slats, or engine pylons for instance) and internal (blade design for turbomachinery) aerodynamics.

Isogeometric analysis 9.2.1 The method

Isoparametric concept consists in the use of the same basis for geometry (i.e. the mesh) and analysis (i.e. the solution space of the numerical method). This concept has been already used in classical finite element analysis where the basis used to approximate the aerodynamic solution field is then used to approximate the given geometry. Therefore, the finite element mesh is only an approximation of the CAD-provided geometry and the resulting approximation can create errors in the solution [3]. Isogeometric analysis is based on the opposite approach. Indeed, a CAD-compatible basis is selected for geometry. It is then used as a basis for the aerodynamic fields.

NURBS as basis functions for geometry

Geometry is constructed from a linear combination of NURBS basis functions. A complete description can be found in [6,3,2]. NURBS functions N p i are described as:

N p i (ξ) = B p i (ξ) • ω i n j=1 B p j (ξ) • ω j ∀ξ ∈ Ξ (9.1)
where:

• Ξ = {ξ 1 , ξ 2 , ..., ξ n+p+1 } is the knot vector, p the polynomial order and n the number of control points (or again the number of degrees-of-freedom) of the geometry parametrization.

• the n B-spline functions of order p (B p i ) are recursively given by the Cox-de Boor formula [7,8,3]:

if p = 0 : B 0 i (ξ) = 1 if ξ i ≤ ξ ≤ ξ i+1 0 otherwise if p ≥ 1 : B p i (ξ) = ξ -ξ i ξ i+p -ξ i B p-1 i (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 B p-1 i+1 (ξ) i = 1, ..., n -1 (9.2)
• ω i (i = 1...n) are weights.

An example of a B-spline basis B p i is shown in Fig. 9.1a. They are polynomials of degree p and are defined on the compact support [ξ i ; ξ i+p+1 ]. The functions B p i can be differentiated p -m times at the knots repeated m times. On the interval [ξ i ; ξ i+1 ], there are p + 1 B-spline functions which are not zero. The NURBS curves C and surfaces S can be defined as follows:

C(ξ) = n i=1 N p i (ξ)P i ξ ∈ Ξ S(ξ, η) = n i i=1 n j j=1 N p i p j ij (ξ, η)P ij (ξ, η) ∈ Ξ × H (9.3)
where N

p i p j ij
are simply given by:

N p i p j ij (ξ, η) = B p i i (ξ) • B p j j (η) • ω ij n i k=1 n j l=1 B p i k (ξ) • B p j l (η) • ω kl (9.4)
An example of NURBS curve is shown in Fig. 9.1b. Note that the curve is interpolatory only at the first and last points. Indeed, B-splines are interpolatory at knots ξ i if and only if their multiplicity m j ≥ p. The fit between the coordinates in the parametric space ξ = (ξ, η) and the coordinates in the physical space X = (x, y) (i.e. the CAD mesh) is given in 2D by:

X(ξ, η) = n i i=1 n j j=1 N p i p j ij (ξ, η)X ij (9.5)
where X ij are the control points defining the CAD mesh (geometry). An illustration of Eq. (9.5) is provided in Fig. 9.2.

Compared to B-spline, the use of NURBS functions is necessary to represent some curve or surface types as conical geometries (circles, ellipses, hyperboles, cones), spheres or revolution surfaces. CAD provides a set a control points with corresponding weights ω i . This allows to represent complex geometric configurations with curved boundaries. In the next section, we show how the CAD-based control points are used to approximate the aerodynamic solution. Figure 9.2: Fit between the parametric space (left) and the physical space (right). (X ij ; Y ij ) are the coordinates of the 5 × 5 control points (red squares). 2 nd order B-splines. Corresponding knot vector is {0, 0, 0, 1, 2, 3, 3, 3} × {0, 0, 0, 1, 2, 3, 3, 3}. Reprinted from [1].

NURBS as basis functions for the aerodynamic flowfield

Following the isoparametric concept, the B-spline/NURBS basis functions are also used for the solution space of the aerodynamic flowfield. The key-point of the isogeometric analysis is to use the CADbased control points as degrees of freedom for the aerodynamic field. Let W (ξ, η) be the conservative aerodynamic field defined on the parametric space Ξ × H (in 2D). The control points related to the aerodynamic field W ij are added to the CAD-based control points X ij (see Eq. (9.5)). This results in the control points P ij = (X ij ; W ij ) and Eq. (9.5) becomes:

X(ξ, η) W(ξ, η) = n i i=1 n j j=1 N p i p j ij (ξ, η) X ij W ij (9.6)
The control points X ij are CAD-provided whereas W ij are the degrees of freedom of the aerodynamics problem to be solved. Nevertheless, both X ij and W ij are at the same location P ij . Note that because of the non-interpolatory nature of the basis,

W (X ij , Y ij ) = W ij .
Regarding 2D Euler compressible system, the stationary equations are solved by a pseudo-time integration:

∂W ∂t + ∂F(W) ∂x + ∂G(W) ∂y = 0 with t -→ ∞ (9.7)
with F and G the fluxes. Following a finite-element formulation, Eq. (9.7) is multiplied by a Bspline/NURBS basis functions N

p i p j ij
and integration is performed all over the computational domain Ω, resulting in the implicit in time weak formulation:

Ω N p i p j ij (x, y) W n+1 (x, y) -W n (x, y) ∆t dΩ = - Ω N p i p j ij (x, y) ∂F n+1 (W) ∂x + ∂G n+1 (W)
∂y dΩ ∀i, j (9.8) After a classical integration by parts and the projection of W on the basis N p i p j ij (see [1] for details), the linear system is inverted with a GMRES algorithm [9]. Note that a stabilizing term (SUPG) for the finite element formulation has to be used to ensure stability of the method [10,11].

Assessment and comparison with classical finite volume methods

The 2D subsonic flow in a nozzle was studied. CAD geometry is shown in Fig. 9.3a with 51 × 11 control points and 3 rd order B-splines. Figure 9.3c represents the Mach number field in the nozzle.

A comparison is performed between the classical 2 nd order finite-volume code (elsA) for two different grids (left: 257 × 65, middle: 51 × 11 points) and isogeometric analysis (right: 51 × 11 control points). The solution computed with the isogeometric analysis-based code (Fig. 9.3c, right) is comparable to those given by the reference solution. On the other hand, the main features of the flow like symmetry with respect to the nozzle axis are missing with elsA 51×11 (Fig. 9.3c, middle), indicating an inaccurate space resolution. Therefore, isogeometric analysis appears to be highly accurate with 3 rd B-splines as basis functions. Only 51×11 degrees of freedom (control points) are needed to reach the same accuracy as a classical finite-volume based code with 257×65 points. Similar conclusions can be drawn regarding the stagnature pressure losses (Fig. 9.3b). However, nuances have to be added to these conclusions. It appears that for the same number of degrees of freedom, isogeometric analysis seems to be more accurate than the classical finite volume method. However, the order p for the B-splines or NURBS functions has to be increased. This has consequences on both the computation time and the stability. Indeed, as introduced in Sec. 9.2.1, the support of the B-spline B p i is the interval [ξ i ; ξ i+p+1 ] whose size increases with p. The direct consequence is the built of a Jacobian matrix which is all the more full as p increases. The cost of the Jacobian matrix inversion, necessary for the resolution of the time implicit linear system, can become prohibitive when the order of the basic functions increases. Another difficulty consists in the adjustment of the SUPG stabilizing term for the high values of p. A difficulty which sometimes proved to be insurmountable.

Conclusions and perspectives

Isogeometric analysis is part of high-performance and high order computing methods. However, its extension to Navier-Stokes equations and its applicability to stiff problems are not obvious. Moreover, the method is only of real interest if it is coupled with a high order grid generator able to take CAD into account accurately. My transfer from DSNA to DMAE and the associated change of topic led me to leave my activities on isogeometric analysis aside. Moreover, at DSNA, emphasis has been placed on other high order methods like Discontinuous Galerkin methods for instance. constraints. Pironneau [10] was the first to use the adjoint equations for design. In the framework of aeronautical CFD, Jameson developed optimal design methods based on optimal control theory. He developed the adjoint formulation for potential flows and Euler equations ( [11,12,13,14]) as well for the Navier-Stokes equation [15]. Complex geometries with complete aircraft configurations dealing with a multiblock implementations were considered [16,17,18,19]. Since then, other people have addressed the problem with the development of adjoint CFD codes for design optimization. A review of the dedicated methods is proposed in [20].

Adjoint methods may also be used as a means for error control in numerical simulations. The adjoint solution may then be seen as an indicator of the sensitivity of the objective function to the local truncation errors in the numerical discretization. Among the methods dedicated to optimal error control are the so-called "goal oriented" mesh adaptation strategies. They were introduced to get satisfactory values of the objective functions at an acceptable cost, using local node displacement and insertion of new points rather than mesh refinement over the entire computational domain. In the context of finite element analysis, important contributions include the work of Johnson et al. [21,22,23], Giles et al. [24], Prudhomme et al. [25], Larson et al. [26], Machiels et al. [27], Hartmann et al. [28,29,30] and Alauzet et al. [31]. In the framework of finite-difference/finite-volume methods, contributions are less numerous.

Basics of the adjoint method

A complete presentation of the adjoint method can be found in [9]. Two formulations are possible for the adjoint method: the discrete one and the continuous one. They are summarized in Fig. 10 both methods, the starting point is the linearization of the original fluid dynamic equations, either written in the continuous form for the continuous adjoint formulation or in the discretized form for the discrete adjoint formulation. For the discrete adjoint formulation, the final step consists in writing the discrete adjoint problem directly from the linearized discrete equations. Regarding the continuous adjoint formulation, two steps are necessary. First, the continuous adjoint problem is formulated from the continuous linearized equations. It is then discretized in a second step. The theoretical basis of the adjoint method is now presented in the context of linear algebra. Indeed, this is the core of the discrete adjoint formulation where the algebraic equations derive from the linearization of the discretized fluid dynamic equations. Moreover, emphasis is put on the discrete adjoint method since this is the one used in the ONERA elsA code optimization solver. The extension to the continuous adjoint formulation is not presented here but can be found in [9] or in [14] where the continuous adjoint method is applied to the Euler equations.

Suppose you compute the aerodynamic flowfield for a given geometric configuration parameterized by several parameters α i and that you want to find the set of α i which optimize (minimize, most of the time) an objective function J . Without loss of generality and to simplify the notations, J will be supposed to be the function of a single parameter α. A change in α may result in a change in the mesh node coordinates (denoted X) and a change in the aerodynamic flow field solution (denoted W ). If the objective function J is written:

J (α) = J (W (α) , X (α)) (10.1)
the gradient dJ /dα is therefore given by:

dJ dα = ∂J ∂W • dW dα + ∂J ∂X • dX dα (10.2)
where X and W are noted in column whereas the gradients of scalars ∂J/∂W and ∂J/∂X are noted in line. The following constraint has to be satisfied simultaneously:

R (W (α), X(α)) = 0 (10.3)
which stands for the non-linear flow field equations (the Navier-Stokes equations for instance). In Eq. ( 10.2), the derivative terms ∂J/∂W and ∂J/∂X are inputs and computed once and for all at the beginning of the optimization process. The term dX/dα is simply calculated by finite differences (f d):

dX dα ≈ dX dα (f d) = X (α + dα) -X (α) dα (10.4)
On the other hand, the term dW/dα is obtained by solving the following linear system derived from the linearization of Eq. ( 10. The right-hand side of Eq. (10.5) is approximated by finite differences (f d):

∂R ∂X dX dα ≈ ∂R ∂X dX dα (f d) = R (X (α + dα) , W (α)) -R (X (α) , W (α)) dα (10.6)
which is a straightforward computation. Indeed, the most expensive operation, namely the calculation of the aerodynamic flow field W , is not influenced by a variation dα of α in Eq. (10.6). Suppose you have n α parameters α. Equation (10.5) has to be solved as many times as n α . Therefore, if the number of parameters n α is quite larger than the number n J of objective functions J to be evaluated, the computation of the gradient dJ /dα with Eq. (10.2) (the so-called direct method) is not efficient. In practice, n α is very large (several hundreds of design variables). The adjoint method allows a computation of dJ /dα whose complexity is independent of n α . It only depends on n J which remains in practice moderate (n J n α ). To do this, the following adjoint equation is solved:

∂R ∂W T λ = - ∂J ∂W T (10.7) 142CHAPTER 10
. ADJOINT METHOD FOR OPTIMAL AERODYNAMIC DESIGN AND ERROR CONTROL which is a linear system whose unknown is the adjoint vector λ (written in column) and which does not depend on n α but on n J . From Eq. (10.2), the computation of the problematic term ∂J ∂W • dW dα can be rewritten using Eq. (10.7) and Eq. (10.5):

∂J ∂W • dW dα = -λ T ∂R ∂W dW dα = λ T ∂R ∂X dX dα (10.8)
Finally, the gradient dJ /dα can be written:

dJ dα = λ T ∂R ∂X dX dα (f d) + ∂J ∂X dX dα (f d) (10.9)
From a practical point of view, Eq. (10.7) is solved by an iterative quasi-Newton method: .

∂R ∂W T (AP P ) λ (l+1) -λ (l) = - ∂R ∂W T (EXA) λ (l) + ∂J ∂W T ( 10

Application of the adjoint method to optimal aerodynamic design

The application of the adjoint method to optimal aerodynamic design was the opportunity for me to conduct technical developments in the optimization solver of the ONERA elsA code. A painful but essential activity is the computation of the (quasi-)exact Jacobian matrix ∂R ∂W T (EXA)

(see Eq. (10.10)) for all the numerical fluxes proposed by the elsA code. And there can be many of them. Unlike the classical Jacobian matrices computed for the implicit time resolution of the Navier-Stokes equations ( ∂R ∂W T (AP P ) for example, Eq. (10.10)) where only the first order fluxes are derived, the derivation of the second order fluxes is necessary for the adjoint method. Any change in the formulation of a numerical flux involves the rewriting of the corresponding derivative terms in the Jacobian matrix. This process is affordable for the inviscid Euler fluxes. However, this is less the case for the viscous fluxes in the Navier-Stokes equations and some symplifying assumptions like the "frozen µ t " hypothesis for the turbulent flows can be done. The "frozen µ t " approximation consists, for the computation of the adjoint vector, in deriving only the first five RANS equations by ignoring the dependence of µ and µ t to the mean field and to the turbulent variable. This approximation brings simplifications but may induce errors on the computation of the gradient vector. To illustrate this, the full linearization of the Spalart Allmaras turbulent model is compared with the "frozen µ t " hypothesis for the computation of the gradient of several objective functions: C Lp (pressure lift coefficient), C Dp (pressure drag coefficient), C Dw (wave drag coefficient), C Df f (far field drag coefficient) and C Dvp (viscous pressure drag coefficient) [6]. A RAE2822 profile is used and the parameter α is the angle of attack. The results are summarized in Tab. 10.1. The use of the "frozen µ t " assumption leads to 10 to 30% errors on the gradient vector computation. The reference is given by the finite difference (FD) computation of the gradients. Reprinted from [6].

Still in the context of numerical methods for the adjoint equation, we proposed a multigrid resolution of the linear system (10.10) similar to the one presented in [32]. The multigrid V cycle between the fine grid h and the coarse one H is described in [33]. Table 4.3 from [33] compares the performances of different methods for the resolution of Eq. (10.10). The test case is the AS28G wing with an inviscid flow. Four different objective functions are studied. The multi-grid resolution significantly reduces the number of quasi-Newton steps needed. This reduction can reach 50% for a three-grid V-cycle. The multi-grid method reduces the CPU time as well.

Regarding the adjoint method implementation in an existing CFD code, the need to ensure the complete derivation of the numerical fluxes (including boundary conditions), is one of the drawbacks of the method. An alternative to "hand derivation" is the use of AD (Automatic Differentiation) software such as TAPENADE 1 to generate the Fortran/C code automatically for the computation of ∂R ∂W T (EXA) (in fact, rather the matrix vector product ∂R ∂W T (EXA) λ (l) from Eq. (10.10)). This option is being used increasingly and we tested it in [33]. However, the structure of the CFD code in which the adjoint method has to be implemented must be thought through since its genesis for an easy integration of the automatic differentiation tools.

Application of the adjoint method to optimal error control 10.4.1 State of the art

At the end of the 1990s, Pierce and Giles introduced adjoint-based error estimation for functions in a very broad framework [34,35]. Consider the linear differential equation Lw = f and suppose the goal to be evaluated is given by (g, w) where (., .) denotes an integral product on a well-posed space Ω (see [1] for more details). The corresponding adjoint problem is given by L λ = g. The error in the estimation of the common goal (g, w) by the approximation solution w h (h denoting the average mesh size) is:

(

g, w) -(g, w h ) = (g, w -w h ) = (L λ, (w -w h )) = (λ, L(w -w h )) = (λ, f -Lw h ) (10.11)
If the adjoint problem itself has also been approximately solved, the error can be expressed as: Pierce and Giles demonstrate that the main error in the previous expression is given by the term (λ h , f -Lw h ) and that in common cases, the order of (λ -λ h , f -Lw h ) is twice the order of the first term. It is noticeable that the error term (λ h , f -Lw h ) in the function of interest is expressed as a weighted sum of the local residual errors of the direct problem with the adjoint variables as the weighting functions. In a series of three articles [36,37,38], Venditti and Darmofal proposed similar formulas for the specific case of finite differences/finite-volume. They presented applications to compressible flow. A Taylor's expansion of the functional output of interest J h about the interpolated coarse-grid solution yields:

(g, w) -(g, w h ) = (λ h , f -Lw h ) + (λ -λ h , f -Lw h ) ( 10
J h (W h , X h ) = J h (W H h , X h ) + ∂J ∂W W H h (W h -W H h ) + O ||W h -W H h || 2 (10.13)
where W is the flow field, X is the volume mesh and R is the residual of the scheme. Two grids are involved: a coarse one H and a fine h. The notation W H h is the coarse-grid flow-field vector reconstructed on the fine grid via some consistent projection operation. Using Eq. (10.7) in this particular context leads to:

λ h | W H h T ∂R h ∂W h W H h = - ∂J h ∂W h W H h (10.14)
and the previous equation can be rewritten:

J h (W h , X h ) = J h (W H h , X h ) -λ h | W H h T ∂R h ∂W h W H h (W h -W H h ) + O ||W h -W H h || 2 = J h (W H h , X h ) + λ h | W H h T R h (W H h ) + O ||W h -W H h || 2 (10.15)
Besides, if the flow computation is not affordable on the fine grid, neither is the solution of Eq. (10.14) for λ h | W H h . The alternative is to replace this adjoint field by the interpolated coarse-grid adjoint:

J h (W h , X h ) ≈ J h (W H h , X h ) + λ H h T R h (W H h ) + λ h | W H h T -λ H h T R h (W H h ) (10.16)
where λ H h represents the coarse-grid adjoint vector reconstructed on the fine grid via some consistent projection operation. The authors recommend to take J

h (W H h , X h ) + λ H h T R h (W H h )
as the function estimate and adapt the mesh by reducing uniformly the error term (last term in the previous equation). See [1] for a description of the different applications of the previous equation.

Later on, Dwight proposed a very different adjoint-based method attached to the scheme of Jameson et al. In a series of two articles [39,40], he considered classical test cases for Euler flows. Computations were run on a hierarchy of grids and for different values of the artificial dissipation coefficients (k 2 , k 4 ). The error for the functions of interest appeared to be mainly due to artificial dissipation. On this basis, the following measure for the approximation error in the scheme of

Goal oriented mesh adaptation based on dJ/dX

From the first use of the dJ/dX method . . . We proposed a method [2,3,4,5,1] based on the total derivative of the goal function J w.r.t. volume mesh coordinates X, denoted dJ/dX. The aim of our approach is to define a new method for finitevolume goal oriented mesh adaptation, in which only one base grid is used (contrary to the approach from Venditti and Darmofal where two grids h and H are used) and which is not restricted to a specific scheme (like for the method from Dwight attached to the Jameson et al. scheme). The gradient dJ/dX is obtained from (see [1] for more details):

dJ dX = λ T ∂R ∂X + ∂J ∂X (10.18)
The computation of dJ /dα simply results in the product of dJ dX (computed by a fast low memory computer since the storage of the volume mesh sensitivities dX/dα is not required) by dX/dα. The vector field (dJ/dX) indicates the variation of the output with the coordinates of the mesh nodes. Nevertheless, the dJ/dX-field may include components orthogonal to the solid walls, which are obviously not usable in the framework of mesh adaptation. In the common case where J is a line integral (in 2D) or a surface integral (in 3D), its total derivative w.r.t. mesh nodes may also include components orthogonal to the integral support which cannot be taken into account during mesh adaptation. This leads to the definition of a projected field denoted P(dJ/dX).

The total derivative P(dJ/dX) can be used to construct global indicators to evaluate the overall quality of a grid for the calculation of the function J. The first indicator (denoted µ J ) is the average of P(dJ/dX) computed with all the nodes of the grid. If large vectors P(dJ/dX) with the same direction are encountered in a zone where the grid is fine, it is not possible to determine an actual nodes displacement that would significantly affect the value of J. This is why another approach consists in multiplying the standard P(dJ/dX) of each node by half the distance to the nearest nearby node. The second indicator, denoted θ J , is the average of this field.

In the framework of goal oriented mesh adaptation based on dJ/dX, we proposed two approaches:

• A node addition method consisting in adding nodes in areas of high ||P(dJ/dX)||-values. This method is suitable for example if the function J cannot be simply correlated with numerical errors.

• A node displacement method well-suited for specific flows for which a functional output is monotonically affected by numerical dissipation. A descent algorithm directly based on P(dJ/dX) or associated with a grid parametrization seeks grids providing better values of the output of interest.

The two methods were applied for both 2D and 3D Eulerian flows computation in [2,3,4,5,1]. At first, only structured O-type grids were used. The extension to unstructured grids was initiated in the PhD study of S. Bourasseau [8]. Regarding the first method (by node addition), a heuristic meshadaptation method, consisting in adding mesh-lines in the zones of large ||P(dJ/dX)|| was successfully applied. The way the nodes (the lines) are added is described in [1]. Basically and if we focus on a 2D configuration, the lines inserted in the coarse grid are interpolated in a reference fine grid. Thereby, a coarse N i × N j grid may be fully defined by the position of its lines (φ i , φ j ), with respect to those of the fine grid (here a 2049 × 2049 grid in 2D). The position of the nodes is then evaluated by a bi-linear interpolation operator. A three-step method is used to add mesh lines to a current mesh. 2D 146CHAPTER 10. ADJOINT METHOD FOR OPTIMAL AERODYNAMIC DESIGN AND ERROR CONTROL Euler subcritical and transonic flows were tested. It was concluded that, contrary to the µ J indicator, the criterion θ J decreases as the grid size is increased and is lower on the stretched grids than on the quasi-uniform grids [1]. This conclusion turned out to be of great importance for the PhD study of M. Nguyen-Dinh [7] since one of the objectives was to define new indicators based on ||P(dJ/dX)||. Similar 3D mesh adaptation by mesh-plane addition were also achieved [1].

Regarding the second method (by node displacement), the visualization of the field P(dJ/dX) gives an insight in the goal-oriented mesh-adaptation issue. If the P(dJ/dX) vector field exhibits a zone of vectors of large magnitude pointing approximately in the same direction, the function of interest J is sensitive to a displacement in this zone. Moreover if the grid is coarse and could be significantly displaced in this area (see Fig. 10.2) then a neighboring acceptable grid would lead to a significantly different value of J and, obviously, a local refinement is needed for a more stable estimation of J. Mesh adaptation without grid parametrization or smoothing was attempted as a first intention. In this case, if steepest descent iterations are used to decrease J, the algorithm simply reads X (l) = X (l-1) -s l P (dJ/dX) (l-1) where s l is a parameter to be adapted. Unfortunately, even for simple configurations, very irregular meshes are obtained after a few iterations of the descent algorithm (Fig. 10.3a), while only a disappointing value of the function of interest is reached. An extended presentation of these results can be found in [3]. After these tests, the projected gradient field P(dJ/dX) has been no longer used directly. Instead, it was combined with suitable parametrizations. To do this, the coarse grid to be optimized is now described by a smooth mapping function (Bézier curves and surfaces) associated with the body-fitted coordinates of the 2049 × 2049 grid. More regular optimized adapted grids are obtained (Fig. 10.3b). Regarding the values of the indicators µ J and θ J for the objective functions C Dp and P a for both the adapted grids and the corresponding quasi-uniform grids, conclusions are similar to those obtained for the line addition method. . . . to more recent developments of the dJ/dX method.

The main conclusions driven from the study in [1] is that the total derivative dJ/dX is a useful quantity for mesh adaptation. However some weaknesses of the presented method were pointed out:

• The proposed indicators (µ J and θ J ) are global and not necessarily the most suitable for the optimization of the objective functions J. That is why the development of more reliable indicators was proposed in the PhD study of M. Nguyen-Dinh [7] and in [41].

• The dJ/dX method cannot be used alone and a grid parametrization is necessary for the mesh displacement method (the use of Bézier curves and surfaces for instance).

• Remeshing strategy has an important impact on the method efficiency.

From these remarks, M. Nguyen-Dinh developed in the second part of his PhD study [7] a more efficient mesh adaptation method for structured grids based on more local refinement indicators. For instance, neighboring points are taken into account. I had no direct involvement in these works due to my transfer from DSNA to DMPE. In his PhD study, S. Bourasseau [8,42] extended the computation of dJ/dX to unstructured and possibly to hybrid grids.

Conclusions and perspectives

If one tries to link this section to the work on icing presented in Sec. I, sensors such as dJ/dα may prove to be effective for assessing accuracy of the computed ice shapes with respect to input data. Indeed, quantities (denoted here J) such as heat transfer coefficient, amount of runback liquid water, sticking efficiency of impinging droplets or solid ice density are very sensitive to parameters (denoted here α) like wall roughness, boundary layer laminar to turbulent transition, local humidity or impinging particle characteristics (melting ratio, velocity, temperature, . . . ). Computation of dJ/dα would then allow to propose an error bar around predicted ice shapes. Error propagation from inputs (those measured by experimenters) towards the final results obtained by the application of the different models could be quantified.

Chapter 11

Conclusions and perspectives

Among my PhD works on the direct numerical simulations of two-phase flows with interface, those on optimization by gradient computation using the adjoint method and those based on methods and models for icing, the common thread for this manuscript did not appear spontaneously. However, whatever my research topic, my activities have always been focused on the two following inseparable points: development of macroscopic models for complex flows and their integration in CFD codes.

The problems encountered in fluid mechanics and energetics are multi-scale and highly non-linear in nature. Thus, the modeling of these phenomena, based on physical considerations and mathematical analysis, may be challenging especially if both small and large scale phenomena are to be modeled at the same time. That is why the great part of my work has been to propose macroscopic models which are supposed to take into account the effects of small scales and their interactions with the larger ones. This approach is driven by the needs of industrial partners for whom the requirements of a DNS computation where all scales are solved are not possible. The derivation of a macroscopic model is conducted in two steps. Firstly, databases, whether experimental or computer-based (DNS), are necessary. Then, and this is the main concern of the approach, the challenge for the modeller is to understand the physical phenomena involved and to extract the relevant parameters to describe them. The space of the parameters in question must be of reasonable size to be easily integrated into a computational code.

In addition to the conclusions and perspectives presented at the end of each chapter of this manuscript, here is a non-exhaustive list of more general perspectives on possible future modeling activities:

Interaction between turbulent and interfacial scales Addressing small interfacial scales remains a challenging task. Indeed, it is at these scales that all the difficulties are combined: low spatial resolution with high curvatures. Often, the front capturing methods are defective at these scales. Although computing capabilities have increased over the last few years, direct numerical simulation with the resolution of all the scales are not yet within the scope of industrial standards. This is why LES modeling of small turbulence scales, especially in the vicinity of deformable interfaces, has a promising outlook. This activity has been continued by S. Vincent after my PhD works. Having left this topic to focus on other projects, I hope to come back to it and renew collaborations on this subject.

SLD and ice crystal icing

Due to the emergence of new certification rules for SLDs and ice crystals for instance, projects have been funded to improve the understanding of the associated physical phenomena. Despite one European project (EXTICE) and two DGAC conventions (PHYSICE and PHYSICE2), no universal model for the sticking efficiency for SLD has been derived. Regarding SLD, two kinds of models exist. The first have a strong theoretical background and the different stages of splashing (from the droplet spreading to the growth of satellite ligaments with corona expansion) are described. Although they are fundamental to the understanding of fine-scale physical mechanisms, unfortunately, these models most often describe single droplet impacts at much lower velocities than those encountered in our aeronautical configurations (∼ 80 m.s -1 ). Furthermore, airflow around the droplet is rarely taken into account. The second kind of models are correlative. Supposedly more representative of reality, they often only fit the database from which they were derived. The reason for this is twofold. Firstly, we must acknowledge that we still have not understood physics of the SLD impacts at high velocity. Secondly, the experimental conditions for which databases are produced are challenging to be generated, measured and calibrated. For instance, large amounts of water (up to LW C = 1 g.m -3 for a velocity of 100 m.s -1 ) must be absorbed by the wind tunnel. Moreover, cloud droplet size distribution and in particular its calibration are tough for icing wind tunnels. Likewise for humidity regulation. Controlling droplet temperature at the impact is complicated, especially for SLDs where thermal equilibrium with supercooled droplets must be ensured. SLDs must be accelerated gradually otherwise they may breakup as a result of air shear. For all these reasons, only a few icing wind tunnels in the world are able to provide experimental conditions well enough controlled to be usable to derive models. But calibration is not everything. Sophisticated measuring devices must be deployed to visualize the relevant physical phenomena. High-speed video are necessary to observe the first moments of the impact. Regarding the amount of water remaining on the wall after a droplet impact at high velocity, the experimental procedure itself is continuously improved. This data is a key factor for sticking efficiency since it increases droplet collection. Another difficulty is water film thickness measurement (in the order of a few tens of microns) on the wall. Moreover, films can be fed by impacting droplets and disturbed by a shearing air flow. Wall roughness characterization and its influence on droplet impact are still being studied. For all these reasons, impact modeling for SLD can only progress through a synergy between experiments focused on the needs of the modeler and numerical simulations that test in real time the derived models on available databases to potentially adjust experiments in return. Regarding the use of DNS, the approach is not mature enough nowadays to derive macroscopic models that could be attractive to industrial partners. However, its relevance is rather in the understanding of small-scale physical phenomena at early stages of splashing and its complementarity with experiments so as to build databases.

Regarding ice crystal icing, the European project HAIC was the first to address the issue. Progress were made on trajectography (drag models and melting rates for flying ice particles), wall impact and accretion. First models, essentially of a correlative nature for those proposed by ONERA, were proposed. All the experimental and modeling difficulties outlined previously for SLDs are obviously relevant for ice crystal icing. To these, are added the difficulty to measure ice crystal shapes and melting ratios. As a first step, the proposed models are most often simplified as for impact where the influence of liquid water on the wall or the nature of the substrate are not taken into account accurately. Although the global mechanism of accretion under ice crystal icing conditions is well understood, particularly inside engines, one of the major gaps at the end of HAIC is the misunderstanding of the role played by liquid water on the wall. This results, for instance, in the right-hand side of the "plateau" (high melting ratios) not being accurately modeled. The water on the wall remains trapped inside the porous ice layer leading to a kind of slushy ice with high erosion rates. Due to a lack of experimental data, the amount of water trapped in the interstices of the porous ice is not available. At most, T. Currie is able to measure a liquid water volume fraction at the wall. Large-scale erosion mechanisms have been identified, namely the increase in the erosion rate with crystal velocity and with the presence of liquid water. On the other hand, fine mechanisms such as discriminating between ice that is eroded and ice that has not sticked due to a reduced sticking efficiency are far from being understood. Within HAIC, the main focus was on unheated walls. The treatment of heated walls which can lead ice crystals to melt when they hit the wall is an additional difficulty that will have to be addressed in MUSIC-haic. The latter has just begun. The objectives are twofold: to improve the models proposed in HAIC and most importantly to deploy them in 3D icing suites.

Ice protection systems and runback modeling

We have just entered a new phase for both models and tools development. Indeed, they are becoming ever more accurate and efficient for the socalled cold configurations, i.e. without ice protection systems. Regarding the latter, we have focused on electro-thermal technologies (ETIPS) at ONERA. Research works on electro-mechanical ice protection systems (EMIPS) have been initiated through the hiring of a post-doc and the submission of a PhD thesis. Back to ETIPS, models dedicated to anti-icing or de-icing must be able to take into account temperature gradients within the substrate. Moreover, an unsteady approach is needed to address phase changes and ice shedding. This is the case for systems operating in de-icing mode where the activation cycles of heating mats are by definition unsteady processes. De-icing is becoming the target for the manufacturers due to reduced power requirements compared to anti-icing mode for which the system is continuously switched on.

ETIPS related scientific challenges have arisen such as modeling of heated boundary layers with enhanced evaporation rates. Topology of the runback liquid water (continuous film, rivulets or isolated droplet) is also of paramount importance. Indeed, depending on it, wet surfaces are different, which has an influence on the thermodynamic balance. The position of the refreezing zone outside the protection system depends on the film topology and dynamics too. This led to the development of ETIPS dedicated models like in MiLeS2D (triple layer) and CRACK2D (ice damage/fracture mechanics).

Significant improvements have been achieved in the modeling of the liquid film breakup into rivulets. Preliminary theoretical explanations for the instability that leads to this fragmentation have been introduced. Preliminary 3D simulations of the development of rivulets have been proposed in the PhD works of J. Lallement. However, there is still room for further improvements like the effect of air shear on film dynamics which has not been fully explored, hysteresis for the contact angle, the consideration of a dirty or rough wall, coupling with heat fluxes leading to film solidification. Moreover, recent runback experiments on heated walls showed that the observed regime (rivulets or isolated drop) is correlated with the wall characteristics (polished aluminum or painted wall). This is why the PhD works of T. Theulier on the runback modeling for isolated drops on a wall have been launched. Since a Lagrangian description of each droplet is not possible, a statistical approach based on the Boltzmann equation will be used.

Icing and classical hard points to be addressed In addition to the specific topics I have just described, there are still some classical hard points to be addressed. Here is a list of some of them. Ice growth process needs local mesh adaptation so that the grid fits the ice shape. Remeshing is acceptable in 2D but not in 3D. To overcome this difficulty, the use of Immersed Boundary Methods (IBM) is being considered. PhD works on the subject are currently in progress at ONERA.

A major difficulty is the computation of boundary layers on rough walls. There are two reasons for this. First of all, the influence of roughness on the boundary layer integral quantities (displacement and momentum thicknesses and shape factor for dynamics), friction coefficient and heat transfer coefficient has to be taken into account. Work has been carried out at ONERA by F. Chedevergne and E. Radenac as part of PHYSICE and PHYSICE2. Secondly, the challenge is to characterize the real rough surface texture by a set of representative parameters such as the equivalent sand grain thickness or the wetted surface for roughnesses. This characterization step is both essential to feed boundary layer models and particularly challenging since it depends on a large number of aero-icing parameters such as ice type (rime or glaze ice) for instance. In addition, experimental characterization of the roughness for real accreted ice remains a tough objective. This gap results in a poor estimation of the heat transfer coefficient, particularly close to the stagnation point. Therefore, heat flux, which is one of the major sink term in the thermodynamic balance, is poorly estimated as a linear function of the heat transfer coefficient. This is particularly acute at temperatures near the freezing point, where such inaccuracy in the heat flux computation can lead to a switch from a glaze to rime ice shape (or vice versa).

A topic is developing: passive ice protection system, not in competition with active ice protection methods but as a complementary one. This activity has to be related to contact angle modeling.

Sensitivity analysis and icing

In the coming years, I would like to share my experience from DSNA on the computation of the sensitivities dJ/dα by the adjoint method and apply it to my icing related problems. Indeed, in the light of what is said above, the evaluation of the sensitivity dh t /dk s of the heat transfer coefficient h t with respect to the equivalent sand grain thickness k s could make it possible to assess the influence of k s on the thermodynamic balance and on the final ice shape. Note that if the number of variables α is not too large, the gradient dJ/dα can be evaluated by direct method without using the adjoint approach.

Icing and numerical tools The development of models has been associated with the writing of an icing suite (IGLOO2D, ETIPS2D, MiLeS2D, CRACK2D, . . . ). The main codes are described in the dedicated chapters of the manuscript. Although it is continuously updated and improved, we can speculate about the future icing suite. As for IGLOO2D, it must be modular. Today, a typical ice shape is obtained from the sequence of stationary solutions for the aerodynamic field, the processing of the dispersed phase, the thermodynamic balance, the computation of the ice shape and the possible re-meshing. A first difficulty, mainly experienced in 3D, consists in making all these codes communicate with each other. Two possibilities are then available. On one hand, file coupling which is simple but costly when reading swap files. On the other hand, shared memory coupling, which is more efficient but requires the use of dedicated couplers. This last approach seems the most promising. For highly unsteady situations (ice shedding for instance) the proposed sequence of stationary solutions may be limiting. Faced to these growing requirements, computational times can increase considerably. To mitigate this phenomenon, so-called fast methodologies are preferred, such as for the computation of the heat transfer coefficient, where an integral formulation for the boundary layer offers a good compromise between efficiency and precision (SIM2D, BLIM2D and BLIM3D).

I have finished with the overview of my research activities since 2009. I hope to maintain this dynamism in the coming years. Icing is a great source of research topics and the near future is promising with the beginning of three icing related PhD thesis at DMPE and two big European research projects MUSIC-haic and ICE GENESIS.

Contents 3 . 1

 31 From supercooled water to ice: basic formulation and application to small droplets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Supercooled Large Droplets (SLD) . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.2 Revisited Model for SLD impact onto a dry solid surface . . . . . . . . . . . . . 24 3.2.3 Direct Numerical Simulations (DNS) of the impact of a single droplet . . . . . 28 3.2.4 Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Figure 3 . 2 :

 32 Figure 3.2: Messinger energy balance. Rime ice configuration (App. C). M ∞ = 0.4, AOA = 1.5 • , P ∞ = 1013.25 hP a, T ∞ = 257.6 K, LW C = 0.31 g.m -3 , M V D = 20 µm, ∆t = 174 s.

  (a) W e = 3255, D = 375 µm, V d = 25 m/s. Time = 0.98 µs. (b) W e = 3255, D = 375 µm, V d = 25 m/s. Time = 4.02 µs. (c) W e = 3255, D = 375 µm, V d = 25 m/s. Time = 20.31 µs. (d) W e = 3255, D = 375 µm, V d = 25 m/s. Time = 36.65 µs. (e) W e = 17500, D = 375 µm, V d = 58 m/s. Time = 0.18 µs. (f) W e = 17500, D = 375 µm, V d = 58 m/s. Time = 3.07 µs. (g) W e = 17500, D = 375 µm, V d = 58 m/s. Time = 6.96 µs. (h) W e = 17500, D = 375 µm, V d = 58 m/s. Time = 12.89 µs. (i) W e = 33300, D = 375 µm, V d = 80 m/s. Time = 0.52 µs. (j) W e = 33300, D = 375 µm, V d = 80 m/s. Time = 2.07 µs. (k) W e = 33300, D = 375 µm, V d = 80 m/s. Time = 3.85 µs. (l) W e = 33300, D = 375 µm, V d = 80 m/s. Time = 6.03 µs.

Figure 3 . 3 :

 33 Figure 3.3: Normal impact of a droplet on a polished aluminum surface. Each line corresponds to a time sequence of an impact for different velocities. Reprinted from [12].

  Definitions of the angle of incidence θ and the impact angle α. (b) against α and |θ|. Comparison between experimental results (symbols) and the present SLD model (solid lines) for different MVDs. MS(1)-317 profile. (c) against α and |θ|. Comparison between experimental results (symbols) and the present SLD model (solid lines) for different MVDs. NACA23012 profile.

Figure 3 . 4 :.Figure 3 . 5 :

 3435 Figure 3.4: against α and |θ|. Comparison between experimental results and the present SLD model for sticking efficiency. Reprinted from [1].

  (a) MS(1)-317 profile at AOA = 0 • . M V D = 137 µm. (b) NACA23012 profile at AOA = 2.5 • . M V D = 154 µm. (c) GLC-305 profile at AOA = 1.5 • . M V D = 92 µm.

Figure 3 . 6 :

 36 Figure 3.6: Collection efficiency β dep against the dimensionless curvilinear abscissa s. Comparison between experiments and numerical results from NASA and ONERA (denoted here "present SLD model"). The "full deposition" (β dep = β imp ) model is shown as well. Reprinted from [1].

  (a) Atmospheric conditions: high Mach number. (b) Altitude conditions.

Figure 3 . 7 :

 37 Figure 3.7: Ice shapes from the DGA database. Reprinted from [1].

Figure 3

 3 Figure 3.8: as a function of θ. The databases of Papadakis, BASE1, BASE2 and BASE3 are represented. Symbols: sticking efficiencies measured experimentally (ONERA) or obtained via Eq. (3.4)(Papadakis). Convex hulls are drawn for small angles of incidence (dashes) and for higher angles of incidence (solid lines). There are two regimes: a regime where tangential effects are predominant and a regime where normal effects are predominant. Reprinted from[12].

Figure 3 . 9 :

 39 Figure 3.9: Water drop impact on a smooth dry solid surface -Yokoï et al. configuration [44]. D = 2.28 mm -water/air/silicon -V = 1m.s -1 . Left: experimental results. Right: DYJEAT 3D Simulation -∆x = 39µm -32 cores -8.3 M cells -T = 35 ms -72 h (SATOR -ONERA). Reprinted from [6].

Figure 3 . 10 :

 310 Figure 3.10: Dynamical contact angle from Yokoï et al.[44] used in the DYJEAT code to compute the case presented in Fig.3.9. Reprinted from[6].

Figure 3 . 11 :

 311 Figure 3.11: Impact at high velocity. M V D = 302 µm, V = 33.6 m.s -1 , water/air, Re = 10147, W e = 4735, Ca = 0.5. Reprinted from [6].

Figure 3 . 12 :

 312 Figure 3.12: Drop impact on a liquid film -Crown formation. Cossali et al. configuration [45]. Reprinted from [6].

Contents 4 . 1 4 . 3 4 . 4 4 . 5 4 . 6

 4143444546 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Ice particle trajectory models . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Geometrical parameters and particle motion model . . . . . . . . . . . . . . . . 4.2.2 Heat exchange and phase change model . . . . . . . . . . . . . . . . . . . . . . Ice particle impingement models . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Impact on unheated walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Impact on heated walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ice particle accretion models . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 The extended Messinger balance . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Erosion model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Model validation: comparison with experimental results . . . . . . . . . . . Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 . 1 :

 41 Figure 4.1: Examples of particles vs. altitude (km) and temperature ( • C) imaged in three size ranges (< 100, 400 -600, > 800 µm by CPI probe on August 22, 1999 [13]. Reprinted from [1].

Figure 4 . 2 :

 42 Figure 4.2: Phase change model for ice crystals. Left: pure ice crystal. Middle: core ice crystal with surrounding liquid water from melted ice. Right: spherical liquid water droplet.

  (a) Experimental and theoretical melting times for all 13 test conditions with error bars -spherical particles only. (b) Experimental and theoretical melting times for all 13 test conditions assuming Φ = 1 -nonspherical particles only. (c) Histogram of the relative differences of the experimental and theoretical melting times for all 13 test conditions assuming Φ = 1 -nonspherical particles only.

Figure 4 . 3 :

 43 Figure 4.3: Comparisons between experimental and theoretical melting times. (a): Spherical particles only. (b) and (c): nonspherical particles only -Φ = 1. Reprinted from [1].

  (a) Experimental and theoretical melting times for 13 test conditions with exemplary error bars. Φ = Φ ⊥ . Nonspherical particles only. (b) Histogram of the relative differences of the experimental and theoretical melting times for 13 test conditions. Φ = Φ ⊥ . Nonspherical particles only.

Figure 4 . 4 :

 44 Figure 4.4: Comparisons between experimental and theoretical melting times. Nonspherical particles only -Φ = Φ ⊥ . Reprinted from [1].

  (a) Evolution of the normal restitution coefficient ξ B n vs. the impact parameter L. Comparison between ONERA fit (solid line)[9], T. Hauk (AGI) experimental results (red circles)[30] and Pan & Render[35] experimental results (blue circles).(b) Evolution of the normal to tangential restitution coefficient ξ F nt vs. the impact parameter L. Comparison between ONERA fit[9] and T. Hauk (AGI) experimental results (largest fragment restitution coefficient)[30].(c) Evolution of the diameter of the biggest fragment vs. the impact parameter L. Comparison between ONERA fit (solid line)[9], T. Hauk (AGI) experimental results (red circles)[30] and Pan & Render[35] experimental results (blue circles).

Figure 4 . 5 :

 45 Figure 4.5: Properties of the secondary (re-emitted) particles. Comparisons between the models from[5, 9] and experimental databases. Reprinted from[30].

Figure 4 . 6 :

 46 Figure 4.6: Crowned cylinder. Stagnation point sticking efficiency variation with η m at M 0 = 0.25 and M 0 = 0.4. P 0 = 34.5kP a (Currie et al. [44]). The function F is plotted with K # = 2.5 (Eq. (4.15)).Reprinted from[4].

  Figure 4.8: "Plateau" effect. Variation of accretion growth with η m for the crowned cylinder at M 0 = 0.25, P 0 = 34.5 kP a and IW C = 6 g.m -3 . t tip is the ice thickness at the stagnation point. See Sec. 4 from [4] for a complete description of the test cases. Reprinted from [44].

Figure 4 .

 4 10b compares the wall liquid mass fraction y l obtained by the model with that measured experimentally by the IPP (y exp l ) for the cylindrical test article at M 0 = 0.25, P 0 = 34.5 kPa and 4 < IW C < 12 g.m -3

  (a) Plateau effect: plateau of almost constant icing severity, which drastically decreases at its left and right limits. (b) Model validation: Plateau effect. Computed ice shapes for run#17, #67, #77, #92 and #101. Comparison between the experimental (red dashed line) and the computed (blue curve) steady state ice thicknesses at the tip.

Figure 4 .

 4 Figure 4.9: "Plateau" effect. See Fig. 4.8 for a comparison with the experimental shapes. Reprinted from [4].

Figure 4 . 10 :

 410 Figure 4.10: Wall liquid mass fraction y l . Reprinted from [4].

  the Cox wind tunnel experiments (runs #9, #10, #19 and #20) performed by Al-Khalil et al. in 2003 [52]: T s ≈ T wb = -7 • C and T s ≈ T wb = -12.5 • C (no melting), M = 0.16, mixed phase, LW C = 0.3 g.m -3 or 0.7 g.m -3 , IW C = 0 g.m -3 , 0.3 g.m -3 or 0.7 g.m -3 , crystal mean mass diameter = 150 -200 µm, • the NASA-NRC wind tunnel experiments (runs #543 and #553) performed by Struk et al. in 2012 [47]: T s = 10.7 • C, T wb = -2 • C and T s = 3.7 • C, T wb = -2 • C (almost no melting due to the negative wet bulb temperature), M = 0.25, mixed phase, LW C = 2.9 g.m -3 , IW C = 7g.m -3 , (a) NACA 0012: Simulation versus experimentside view ice accretion after 120 seconds of icing at T W C = 12 g/m 3 and LW R = 0.28. Variation of temperature. (b) NACA 0012: Simulation versus experimentside view ice accretion after 120 seconds of icing at T W C = 12 g/m 3 and LW R = 0.12. Variation of temperature. (c) Cylinder: Simulation versus experiment -side view ice accretion after 120 seconds of icing at T W C = 12 g/m 3 and LW R = 0.12. Variation of temperature. (d) NACA 0012: Simulation versus experimentside view ice accretion after 120 seconds of icing at T ∞ = 0 • C, T W C = 12g/m 3 for a) and b) and T W C = 6.4 g/m 3 for c). Variation of the liquid water ratio LW R. (e) Cylinder: Time history of leading edge accretion angle ϕ 0 . Variation of the liquid water ratio LW R. (f) NACA 0012: Simulation vs. experiment -Side view ice accretion after 120 seconds of icing at T ∞ = 0 • C and LW R = 0.2. Variation of total water content.

Figure 4 . 11 :

 411 Figure 4.11: Influence of the temperature T ∞ , LW R and T W C on ice shape. Reprinted from [3].

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.2 Overview of IGLOO2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 IGLOO2D assessment (App. C) . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.4 Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 5 . 1 :

 51 Figure 5.1: IGLOO2D organization chart. Reprinted from [1].

  (a) Structured grid. (b) Unstructured grid. (c) Hybrid grid with a layer of structured quadrangles near the profile and unstructured triangles far from the profile. (d) Multi-element airfoils (NHLP).

Figure 5 . 2 :

 52 Figure 5.2: Different types of grids obtained with MESH2D. Reprinted from [1].

5 . 3 .Figure 5 . 3 :

 5353 Figure 5.3: IGLOO2D computational loop from the clean profile to the final iced shape. Reprinted and adapted from [1].

Figure 5 . 4 :

 54 Figure 5.4: Some results of rime ice shape computations with IGLOO2D.

Figure 5 . 5 :

 55 Figure 5.5: Ice shape comparison between LEWICE and IGLOO2D. Experimental results provided by NASA are depicted using solid plain grayscale. The darker the shape, the most likely the experimental ice shape. The experimental ice shapes are obtained from a statistical average (several trials of the same experiment or spanwise average) showing the experimental repeatability. The input conditions are given in Tab. 2 from [1]. Reprinted from [1].

6 . 1 6 . 2 A 6 . 3 A 6 . 4 6 . 5

 6162636465 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 triple layer based software as a substitute for the Messinger approach: MiLeS2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.2.1 Physical modeling and equations . . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.2.2 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.2.3 Assessment and basic computations . . . . . . . . . . . . . . . . . . . . . . . . 89 model for electro-thermal ice protection system . . . . . . . . . . . . . . 90 6.3.1 Thermal modeling: ETIPS2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.3.2 Coupling between MiLeS2D and ETIPS2D with a Schwarz algorithms . . . . . 91 Mechanical modeling and ice shedding: CRACK2D . . . . . . . . . . . . . . 93 Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

  -layer solver solid ice block (white) over a static liquid water film (blue) (a) Coupling between MileS2D and ETIPS2D. Reprinted from [1]. meshing, aerodynamics and trajectory. See Sec. 5 for a complete description (b) IPS solvers in IGLOO2D. MiLeS2D is a more advanced alternative to MESSINGER2D for accretion (see Fig. 5.3).

Figure 6 . 1 :

 61 Figure 6.1: Electro-thermal ice protection systems in IGLOO2D.

  Illustration of the three layers: h f , h i and h s are the heights of the running film, ice and static (melted) film. The a k 's and b k 's represent the positions of the boundaries between each layer in the z direction (normal to the surface).

Figure 6 . 2 :

 62 Figure 6.2: Three-layer model in MiLeS2D. Reprinted from [4].

Figure 6 . 3 :

 63 Figure 6.3: Comparison between x and y components of the temperature gradient at t = 30s. Arrows represent -∇T . The dashed line represents the position of the melting front. Reprinted from [4] (Appendix A).

  (a) A small region of ice is melted. The resulting liquid water runs downstream. (b) Predicted runback ice shapes after 50s (obtained with a grid size of ∆s = 10 -4 m).

Figure 6 . 4 :

 64 Figure 6.4: Computational results obtained with MiLeS2D during unsteady heating of a leading edge. Concept of runback ice build-up (ice ridges). Reprinted from [4].

  (a) Zone distribution example. (b) Layer distribution example.(c) Multi-layered structured grid example.

Figure 6 . 5 :

 65 Figure 6.5: Zone and layer distribution example. Structured grid. Reprinted from [2].

Figure 6 . 6 :

 66 Figure 6.6: Non-overlapping 1D geometrical domains Ω 1 and Ω 2 used for the Schwarz method.

  (a) Non-optimal coupling coefficients ω 1 and ω 2 .(b) Optimal coupling coefficients ω 1 and ω 2 .

Figure 6 . 7 :

 67 Figure 6.7: Speed of convergence of the Schwarz method. Influence of the coefficients ω 1 and ω 2 . u stands for the temperature T . Blue (resp. red): convergence criterion based on the continuity of the temperature field (resp. temperature gradient) at the interface x = 0. Reprinted from [1].

  (a) Electro-thermal ice protection system operating illustration. (b) Geometrical illustration. (c) Schematic pressure distribution.

Figure 6 . 8 :

 68 Figure 6.8: Proposed scenario for ice shedding. Reprinted from [2].

  Clear fractured state through the ice.(c) No fracture if L f /L t is too low.

  (a) Operating cycle of the heating mats (ignition). (b) Ice shape after 100 s. (c) Ice shape after 110.5 s. (d) Ice shape after 160 s.

Figure 6 . 10 :

 610 Figure 6.10: Operating scheme of a complete electro-thermal ice protection system. Reprinted from [1].

Figure 7 . 1 :

 71 Figure 7.1: Thin film configuration. Reprinted from [5].

Figure 7 . 2 :

 72 Figure 7.2: Definition of e f ilm for asymptotic regimes h = 0 and h R.In this figure, for the sake of clarity and without loss of generality, the film is assumed to be in static equilibrium on a horizontal flat plate (e k = 0, Eq. (7.11)). Reprinted from[1].

  (a) Configuration of a steady liquid film on a horizontal flate plate. θ s is the static angle and W s the spreading. h puddle is the film thickness at equilibrium. (b) Influence of h on the global film shape. (c) Influence of h : focus near the stagnation line.

Figure 7 . 3 :

 73 Figure 7.3: Spreading of a liquid film on a horizontal flate plate. Static steady state. Influence of h . The spatial step ∆x is set. Reprinted from [1].

  (a) Film injected at a constant mass flow rate on an inclined flate plate. h N u is the thickness of the Nusselt flow. h bulk is the thickness of the advancing bulk front. (b) Film profiles for different values of b. θ s = 0 • (fully wetting). t = 1 s. (c) Film profiles for different values of b. θ s = 38 • (partially wetting). t = 1 s.

Figure 7 . 4 :

 74 Figure 7.4: Influence of the slip length b (Eq. (7.20c)) on the film profiles. Comparison between fully and partially wetting configurations. The spatial step ∆x and the thickness h are set so that h ∼ h N u /10 and ∆x h . Reprinted from [1].

  (a) Influence of b on h N u (far from the contact line) and h bulk (advancing bulk front). (b) Influence of b on the relative error r b between h bulk and h N u .

Figure 7 . 5 :

 75 Figure 7.5: Fully wetting film. Influence of b on the thicknesses h bulk and h N u . The spatial step ∆x is set. Reprinted from [1].

Figure 7 . 6 :

 76 Figure 7.6: Partially wetting film. Influence of b on the thicknesses h bulk and h N u . The spatial step ∆x is set. Reprinted from [1].

  (a) Numerical simulation. θ s = 68 • . (b) Experimental results. Q = 9.17 × 10 -6 m 3 .s -1 . (c) Numerical simulation. θ s = 80 • . (d) Numerical simulation. θ s = 68 • . (e) Experimental results. Q = 1.25 × 10 -5 m 3 .s -1 . (f) Numerical simulation. θ s = 80 • .

Figure 7 . 7 :

 77 Figure 7.7: Thickness profile at t = 0.3 s of a liquid water film injected on a vertical flate plate as a function of the volume flow rate Q. The static contact angle θ s is measured experimentally and estimated between 60• and 80 • . h is chosen so that h /h bulk ≈ 1/18. The pinching length is symbolized by the white solid line. Reprinted from[1].

Figure 7 . 8 :

 78 Figure 7.8: Numerical simulations based on Johnson's experiments[23]. Re = 0.52. The configurations are summarized in[1]. Reprinted from[1].

  (c) Bell et al. test case[27]. 64 2 grid. Time evolution of the mass of the initial circle for different interface capturing methods. Acronyms are available in[1].

Figure 8 . 1 :

 81 Figure 8.1: Test cases from Zalesak[24] and Bell et al.[27]. The numerical methods used are detailed in[1]. Reprinted from[1].

  (a) Liquid sheet in a HIT flow for four different W e numbers at different times t/T e . (b) Anisotropy tensor a uu and a ww . Anisotropy is enhanced near the interface (z/δ → 0). (c) 3D energy spectra. Solid line: W e = ∞ (no surface tension). Dashed line: W e < ∞. Black: t/T e = 3. Red: t/T e = 6. Blue: t/T e = 13. Green: t/T e = 20. Orange: t/T e = 30. 512 3 grid.η is the Kolmogorov space scale at t/T e = 3. Energy is transferred from large (kη 1) to small scales (kη ≈ 1) when the interface surface tension is not zero. t is time and T e is the Eulerian time scale or eddy turnover time based on the large scales of turbulence.

Figure 8 . 2 :

 82 Figure 8.2: Reprinted from [1].

Figure 8 . 3 :

 83 Figure 8.3: Comparison of the filtered terms τ k and τ Φ with their proposed models τ mod k and τ mod Φ . t is time and T e is the Eulerian time scale or eddy turnover time based on the large scales of turbulence. Reprinted from [1].

9 . 3

 93 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 9.2 Isogeometric analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 9.2.1 The method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 9.2.2 Assessment and comparison with classical finite volume methods . . . . . . . . 135 Conclusions and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

(a) 3

 3 rd order B-spline basis. n = 9. The open knot vector is given by: Ξ = {0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 6, 6, 6}. (b) A NURBS curve in R 2 (solid black). Control points are represented (red squares). The curve is interpolatory only at the first and last control points.

Figure 9 . 1 :

 91 Figure 9.1: (a): B-spline basis ; (b): NURBS curve. Reprinted from [1].

  (a) CAD geometry of the nozzle. 3 rd order Bsplines, 51 × 11 control points (black squares). (b) Stagnation pressure evolution along the down wall of the nozzle. (c) Mach number in the nozzle. Left: classical 2 nd order finite-volume elsA code, 257 × 65 grid. Middle: classical 2 nd order finite-volume elsA code, 51 × 11 grid. Right: isogeometric analysis, 3 rd B-splines, 51 × 11 control points.

Figure 9 . 3 :

 93 Figure 9.3: Nozzle test case. Reprinted from [1].

Figure 10 . 1 :

 101 Figure 10.1: Continuous (in red) and discrete (in blue) adjoint methods. Reprinted from [9].

  .10) so that λ (l) -→ l→∞ λ. The matrix ∂R ∂W T (AP P ) is the approximated Jacobian matrix whereas ∂R ∂W T (EXA)is the (quasi-)exact Jacobian matrix. See next section (Sec. 10.3) for the difference between ∂R ∂W T (AP P )

Figure 10 . 2 :

 102 Figure 10.2: NACA0012, subcritical flow conditions. Examination of -P(dC Dp /dX) (left) and P(dP a /dX) (right) on quasi-uniform 257 × 257 grid. Top: general view. Down: leading edge. Respective scales of arrows are 20 (top figure), and 1 (down). Reprinted from [1].

Figure 10 . 3 :

 103 Figure 10.3: Mesh adaptation. Dashed: initial grid ; Solid: adapted grid for C Dp . Reprinted from [7].

Table 3 .

 3

					ṁL dep • H L dep	Qc	ṁevs • H evs
	ṁL dep	ṁevs
				ṁrb,in • H rb,in	ṁrb,out • H rb,out
	ṁrb,in	ṁacc	ṁrb,out
	(a) Mass balance.			Qheat
					(b) Energy balance.
		Figure 3.1: Messinger balance[8]. Reprinted from [9].
	Mass balance (kg.m -2 .s -1 )		Energy balance (J/kg)
	notation definiton		notation definition
	ṁL dep	liquid water deposed surface	H L dep	liquid water deposed mass en-
		mass rate			thalpy
	ṁevs	evaporated/sublimated sur-	H evs	evaporated/sublimated mass
		face mass rate			enthalpy
	ṁrb,in	upstream runback surface	H rb,in	upstream runback mass en-
		mass rate			thalpy
	ṁrb,out	downstream runback surface	H rb,out	downstream runback mass en-
		mass rate			thalpy
	ṁacc	surface ice growth rate	Qc Qheat	convective heat flux heating flux

1: Terms of the Messinger mass and energy balances (Eq. (3.1)) detailed in Tab. 3.1. Historically introduced by Messinger [8], the mass and energy balances can be written: d dt m L = ṁL dep + ṁrb,inṁrb,outṁevsṁacc d dt (U + E k ) -ṁL dep H L dep -ṁrb,in H rb,in + ṁevs H evs + ṁrb,out H rb,out = Qc + Qheat (3.1)

  3.1, where the superscript L has been omitted and where ṁimp is the impacting liquid water surface mass rate). ṁdep stands for the amount of liquid water remaining on the wall. ṁimpṁdep is the amount of water re-emitted from the wall (by splashing). Regarding the properties of the re-emitted droplets, see Sec. 3.2.2.Papadakis et al.

  .7. The particle deposited mass rate ṁdep derives from the

		Solid ice
	ṁimp	Liquid water
	ṁimp -ṁdep
	ṁevs
	ṁrbi	ṁer
		ṁrbo
	ṁL acc + ṁS

acc Figure 4.7: Messinger mass balance. Extension to ice crystal icing. ṁimp , ṁdep , ṁevs , ṁer , ṁrbi , ṁrbo , ṁL acc and ṁS

Table 10 .

 10 1: Comparison between the "frozen µ t " assumption and the full linearization of the Spalart Allmaras turbulent model. RAE2822 profile.

	function	FD	"frozen µ t " assumption full linearization value error value error
	dC Lp /dα	1.703e-01 1.105e-01	35%	1.687e-01 0.9%
	dC Dp /dα	5.688e-03 4.741e-03	17%	5.703e-03 0.3%
	dC Dw /dα 2.272e-03 1.658e-03	27%	2.251e-03 0.9%
	dC Df f /dα 4.092e-03 3.606e-03	12%	4.069e-03 0.5%
	dC Dvp /dα 1.820e-03 1.948e-03	7%	1.818e-03 0.1%

  Jameson et al. was proposed:

	k 2	dJ dk 2	+ k 4	dJ dk 4	(10.17)

The field of local indicator for dissipation error is used as a mesh refinement indicator and J -k 2 dJ dk 2k 4 dJ dk 4 is considered as the corrected output value.

Rime ice refers to opaque, rather brittle and low density ice obtained by fast freezing at low temperatures. On the other hand, glaze ice is a smooth, transparent and homogeneous ice obtained at higher temperatures.

Supercooled large droplets (SLD), as suggested by their name, are supercooled droplets whose diameter is greater than 50µm. The two DGAC conventions PHYSICE and PHYSICE2 were the opportunity for me to focus on the modeling of SLDs. This was the context for writting one paper[1] and several conference papers or oral presentations[2, 3, 4, 5, 6]. I have been involved in the PhD works of T. Xavier supervised by D. Zuzio and J.L. Estivalezes. The subject is Direct Numerical Simulation (DNS) of the droplet impact on a wall.

Acknowledgments

This part is composed of two chapters:

• Chapter 6: description of the models and tools dedicated to electro-thermal ice protection systems (ETIPS).

• Chapter 7: wettability, partially wetting films and rivulet modeling. The objective of this section is to summarize briefly my PhD works conducted between 2006 and 2009 at ONERA. These works, entitled Développement d'une approche de type LES pour la simulation d'écoulements diphasiques avec interface. Application à l'atomisation primaire (Development of a LES approach for the simulation of two-phase flows with an interface. Application to primary atomization), are presented in [1]. They were the subjects of three papers [2,3,4] and several conference papers [5,6,7,8,9,10,11,12].

Part III

Other works

Introduction

Whereas Large Eddy Simulation (LES) of single-phase flows is already widely used in the CFD world, even for industrial applications, LES of two-phase interfacial flows, i.e. two-phase flows where an interface separates liquid and gas phases, still remains a challenging task. The main issue is the development of subgrid scale models well suited for two-phase interfacial flows. The aim of my PhD works was to generate a detailed DNS database of incompressible two-phase interfacial flows in order to clearly understand interactions between small turbulent scales and the interface separating the two phases. The different subgrid scale terms were derived from an a priori analysis (filtering) of this DNS database. To do this, different numerical techniques were firstly compared to find which methods are well adapted to deal with turbulent configurations where large interfacial deformations occur (Sec. 8.2). Then, the interface/turbulence interaction were studied in the configuration where the interface is widely deformed and where both phases are resolved by DNS (Sec. 8.3). The interaction between an initially plane interface and a freely decaying homogeneous isotropic turbulence (HIT) was studied. The densities and viscosities were the same for both phases in order to focus on the effect of the surface

Chapter 10

Adjoint method for optimal aerodynamic design and error control This chapter is focused on the adjoint method which was one of my research topics during my first years at the DSNA (ex DAAA) department at ONERA. In a first part (Sec. 10.2), the basics of the adjoint method is presented. Especially, the discrete formulation is detailed. In a second part (Sec. 10.3), the application of the adjoint method to optimal aerodynamic design is outlined. This activity was an opportunity for me to conduct numerical developments in the optimization solver of the ONERA elsA code. The third part (Sec. 10.4) is dedicated to the use of the adjoint technique for error control. In particular, we developped a goal oriented mesh adaptation technique using the total derivative of aerodynamic functions with respect to mesh coordinates (the so-called dJ/dX technique). This approach was the subject of one paper [1] and several conference papers and oral presentations [2,3,4,5,6]. It is based on the PhD studies of M. Nguyen-Dinh [7] and S. Bourasseau [8]. I co-advised the first part of their works with J. Peter before my transfer to the DMPE department.

Introduction

The adjoint method is involved in optimal aerodynamic design and optimal error control. Regarding optimal design, adjoint solutions have been used for a long time [9]. They aim at providing the linear sensitivities of an objective function (like lift or drag for instance) to a large number of design variables which parametrise the shape. These sensitivities can be used to drive an optimization procedure which simply relies on the use of descent algorithms (steepest descent for instance). The word "optimal" has to be understood as the way to find the geometry which minimises an objective function given some